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We analyze the equilibrium properties of a weakly interacting, trapped quasi-one-dimensional Bose 
gas at finite temperatures and compare different theoretical approaches. We focus in particular on 
two stochastic theories: a number-conserving Bogoliubov (ncB) approach and a stochastic Gross- 
Pitaevskii equation (sGPe) that have been extensively used in numerical simulations. Equilibrium 
properties like density profiles, correlation functions, and the condensate statistics are compared to 
predictions based upon a number of alternative theories. We find that due to thermal phase fluctu- 
ations, and the corresponding condensate depletion, the ncB approach loses its validity at relatively 
low temperatures. This can be attributed to the change in the Bogoliubov spectrum, as the conden- 
sate gets thermally depleted, and to large fluctuations beyond perturbation theory. Although the 
two stochastic theories are built on different thermodynamic ensembles (ncB: canonical, sGPe: grand- 
canonical), they yield the correct condensate statistics in a large BEG (strong enough particle inter- 
actions). For smaller systems, the sGPe results are prone to anomalously large number fluctuations, 
well-known for the grand-canonical, ideal Bose gas. Based on the comparison of the above theories 
to the modified Popov approach, we propose a simple procedure for approximately extracting the 
Penrose-Onsager condensate from first- and second-order correlation fimctions that is computation- 
ally convenient. This also clarifies the link between condensate and quasi-condensate in the Popov 
theory of low-dimensional systems. 

PACS numbers: 03.75.Hh; 67.85.-d; 67.85.Bc; 05.10.Gg; 42.50.-p 
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I. INTRODUCTION 

Following the first observation of Bose-Einstein con- 
densation of dilute gases, experimental and theoreti- 
cal efforts were mainly focused on the fundamental 
properties of such degenerate quantum gases, includ- 
ing spatial and momentum distributions, and collec- 
tive excitations [1 2J. Mean field theory was initially 
found to be impressively successful in most cases. In- 
deed, at temperatures well below the phase transition, 
nearly all atoms occupy one wave function that satisfies 
a nonlinear Schrodinger equation, the celebrated Gross- 
Pitaevskii equation (GPe). The nonlinear term describes 
the mean field potential experienced by the atoms due 
to two-body interactions. In the language of quantum 
field theory, the GPe yields a zeroth-order approxima- 
tion to the full matter-wave field where both the non- 
condensed component of the gas and quantum fluctua- 
tions are neglected. 

The Bogoliubov theory provides an improved analy- 
sis of this system by including small fluctuations around 
the condensate wave function. Its predictions include, 
e.g., the spectrum of collective excitations, the quantum 
depletion of the condensate, and correlation functions 
at both zero and nonzero temperature. We focus in this 
paper on a one-dimensional trapped gas as a model for 
a weakly interacting quasi-one-dimensional system con- 
fined tightly in the radial direction. In such a system, the 
contribution of low-energy modes is significant. From 
Bogoliubov theory, these modes mainly affect the phase 
of the matter-wave field, the density fluctuations being 
relatively weak. As a consequence, there is no Bose 
condensate in the homogeneous limit. Still, a so-called 
quasi-condensate can be identified where long-range co- 
herence manifests itself in the suppression of density 
fluctuations, while the phase is correlated only over dis- 
tances smaller than the system size fSHSI. The situation 
is similar to a "fragmented" condensate where several 
low-energy modes appear with comparable weight |9|. 

The (quasi-)condensateand atoms in excited states 
("thermal cloud") are often treated as two subsystems 
that are coupled to each other by scattering processes 
that exchange particles and energy. Approaches based 
on such a splitting between condensate and thermal 
d5mamics lead to a generalized GPe for the conden- 
sate dynamics that differs from its T = counter- 
part through the inclusion of the thermal cloud mean 
field (Hartree-Fock potential). In addition, a source 
term may describe the scattering of particles between 
the condensate and thermal cloud llOi 1111 . The ther- 
mal cloud itself is described by a quantum Boltzmann 



equation [1^-14] self-consistently coupled to the con- 
densate 1 15 - 19 1 . In its kinetic formulation, the resulting 
self -consistent, coupled Gross-Pitaevskii-Boltzmann ap- 
proach, which extends earlier work by Kirkpatrick and 
Dorfman |20|, and Eckern [21 1, is often referred to as 
the "ZNG" scheme within the context of trapped atomic 
gases [15, 22 J. This method reproduces the two-fluid 
hydrodynamics in the collisional, hydrod5mamic regime 
| |23ll24| , and has been tested successfully against experi- 
ment for collective modes [25-27J and macroscopic exci- 
tations Il28ll29ll . Since the ZNG approach is numerically 
formulated in a purely 3-dimensional context, we shall 
not be considering it further within the present work, 
based on purely one-dimensional simulations. 

Despite their elegant formulation, kinetic theories 
based on mean field potentials have the drawback in 
lower dimensions of not fully capturing the large phase 
fluctuations in the quasi-condensate. In addition, so- 
called anomalous averages (or pair correlations) in the 
thermal cloud, typically omitted in such approaches, 
are expected to become particularly relevant at lower 
dimensions. It is not entirely clear, however, how to 
obtain a gapless excitation spectrum for the system in 
the homogeneous limit, as required by the Hugenholtz- 
Pines theorem ll6l[7l l30H34| . A modified mean field the- 
ory for low-dimensional quasi-condensates was devel- 
oped by Stoof 's group and one of the present co-authors 
jUJIZl, building on previous path-integral approaches pi- 
oneered by Popov [3 5|. In this "modified Popov the- 
ory", the infrared divergences due to phase fluctuations 
are systematically removed, leading to a gapless, con- 
vergent and computationally convenient scheme that 
applies in all dimensions, for homogeneous and trapped 
systems. This approach has been used in particular to 
construct the phase diagram of weakly-interacting ID 
systems [35 1, and to investigate the interplay of density 
and phase fluctuations \36}. 

The quasi-condensate dynamics at nonzero tempera- 
ture is a challenging problem as the Bogoliubov approx- 
imation becomes invalid, at least at large times [37,,38J , 
and large thermal phase fluctuations have to be taken 
into account even at low temperatures where density 
fluctuations are small [39 J. Classical field methods have 
been developed to simulate numerically those modes of 
the system that feature a macroscopic occupation. These 
methods rely upon the observation that for highly oc- 
cupied modes, the field operator can be replaced by a 
classical complex field which evolves in time according 
to the GPe [^-42 J. This description extends the T = 
GPe, by adding stochastic elements that describe fluc- 
tuations of the (quasi-)condensate modes; these may be 
further coupled to the thermal cloud where the mean 
occupation numbers are small. (In fact, too small to be 
treated in the classical field approximation, and more 
appropriately described by quantum Boltzmann equa- 
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tions t^)-) Within the class of classical field techniques, 
we mention the projected GPe (pGPe) Ii41n44j , the trun- 
cated Wigner (tW) method 1^ I45H47I the stochastic 
Gross-Pitaevskii equation (sGPe), when implemented in 
the classical limit |43, 48-53|, and closely related classi- 
cal field methods [42, 54 J. Hybrid simulation techniques 
were also recently developed that attempt to go beyond 
the classical limit [55 - 58]. 

These stochastic approaches, the relation between 
them and other kinetic theories and their respective ap- 
plications have been reviewed in Refs.f52l l53ll59ll60l . A 
key appeal is that they provide an approximation to the 
full distribution function of the ultracold gas and give 
access to physics beyond the mean field. They have been 
used, e.g., to probe the large critical fluctuations near 
the phase transition ll6Tti63] , to study dynamical effects 
of fluctuations on condensate growth H48l |49l |64ll and 
macroscopic excitations 165- 691 . Another quantity of in- 
terest is the counting statistics of the condensate mode 
Il63l [70H73l , which is analogous to the photon number 
distribution in quantum laser theory [74]. 

As the use of classical field simulations becomes more 
widespread, a quantitative study of their relative prop- 
erties is essential. The main purpose of this paper is 
to initiate such a quantitative study by comparing two 
methods that can be implemented with reasonable ef- 
fort, each of which generates an ensemble of stochastic 
initial states to mimic a finite-temperature Bose gas at 
equilibrium in a trap. For simplicity, we focus on a one- 
dimensional, weakly interacting system. 

Number Conserving Bogoliubov (tWncB): The first 
method is a Bogoliubov approach in which the to- 
tal atom number N is conserved (formulated within 
the canonical ensemble) [46 77-81]. This method 
has been used as a starting point for d5mamical cal- 
culations within the truncated Wigner approximation 
(see Ref. |60| for a review). We henceforth denote 
this by tWncB. The tWncB field contains both conden- 
sate and non-condensate modes, calculated from the 
Bogoliubov-de Gennes equations. The modes ampli- 
tudes are sampled to capture both quantum and ther- 
mal fluctuations. We adopt here a formulation devel- 
oped in a low-temperature expansion around the Gross- 
Pitaevskii mean field 1461 |80l , assuming that Nt^, the 
number of atoms in non-condensate modes, is small 
compared to N. 

Stochastic Gross-Pitaevskii Equation (sGPe): The sec- 
ond method is the stochastic Gross-Pitaevskii equa- 
tion (sGPe), which prepares a grand-canonical ensem- 
ble dynamically by simulating a Langevin equation (see 
Refs.im IMl |49] for details of the scheme used here). 
We consider it here within the classical approximation 
where the mode occupations are large. The stochas- 
tic field in the sGPe represents the low-lying modes of 
the field which are coupled to a thermal cloud, treated 



as a heat bath. The exchange of particles and energy 
through incoherent scattering processes between these 
two sub-systems is represented by a damping term and 
the Langevin force in the sGPe [43 1. 

We wish to show that there exist regimes where the 
two methods are equivalent despite the physically very 
different pictures behind them. To this end, we calculate 
and analyze relevant observables like density profiles, 
spatial correlation functions, and condensate statistics. 
Where feasible, we also compare to alternative finite 
temperature theories, as detailed in Table |l] This study 
is by no means complete (e.g., there is no comparison 
to 'ZNG' or pGPe), but we believe that it represents an 
important step towards benchmarking commonly used 
simulation methods for finite-temperature Bose gases. 
Other comparisons undertaken to date are summarized 
inRef.[82]. 

More specifically, we explain the origin of discrepan- 
cies between the two methods considered here, build- 
ing on previous investigations of the validity conditions 
of the tWncB [38]. We find in particular that the low- 
temperature (or small A^th) expansion behind the tWncB 
breaks down quite early, as the temperature T increases 
towards the characteristic temperature for phase co- 
herence within a trap [4]: 



(1) 



(cj is the trap frequency and ^ the chemical potential) 
where T^. is the critical temperature for Bose-Einstein 
condensation in an ideal trapped Bose gas in ID [76] . 
This failure, that also happens when the theory of 
Ref.[71] is applied to an interacting, trapped gas, is at- 
tributed to a distribution function for TVth that is broad- 
ened by thermal fluctuations beyond the limit set by the 
total number of particles N. In addition, these fluctu- 
ations are overestimated because of spurious contribu- 
tions from phase fluctuations. We have observed that 
the breakdown of the tWncB approach is not completely 
"cured" by propagating the stochastic ensemble of wave 
fields under the GPe. 

At low temperatures and for smaller systems, we 
have found large number fluctuations in the sGPe re- 
sults, in particular in the counting statistics of the con- 
densate. This is related to the anomalous number fluctu- 
ations of the ideal gas in the grand-canonical ensemble 
Il70ll83ll84ll . This feature does not occur for the canonical 
tWncB method, and it is removed in larger condensates 
where particle interactions become important. We em- 
phasize that this agreement illustrates how moments of 
the quantum field of very high order are correctly repro- 
duced by the stochastic approaches. 

In addition, we discuss the influence of the ther- 
mal (non-condensate) density nth(^) and the so-called 
anomalous average (or pair correlation) m{z) of the 
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Condensate statistics 
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Svidzinsky and Scully 
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Pair anomalous average 
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(T — 0) Bogoliubov theory 
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/ 




modified Popov 
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Table I: Summary of 'benchmark' theories used in analyzing the stochastic approaches. A tick in the 'ncB' or 'sGPe' columns 
indicates inclusion in the comparison. [Eq.([l|] is the critical temperature above which the system's coherence length is smaller 
than its size due to phase fluctuations li4J. Tc [Eq.jTJ] is the critical temperature for an ideal Bose gas in a one-dimensional 
harmonic trap l76l . These temperatures are illustrated in Figlllin relation to our parameters. 



non-condensate field, by analyzing their back-action on 
the shape of the (Penrose-Onsager) condensate den- 
sity. This anomalous average is related to both a renor- 
malization of the particle interactions due to the back- 
ground field, and to the Landau and Beliaev damping 
of condensate excitations (together with triple averages) 

E llQl El [13 123 IMl IMl ESMl . 

Following on from this, we discuss the connection 
between the condensate mode, as obtained by ap- 
plying the Penrose-Onsager criterion, and the quasi- 
condensate, as predicted by the modified Popov theory 
of Andersen et al. \6.7\. In stochastic theories, the con- 
densate mode is commonly extracted a posteriori from 
the total matter field ensemble, which can become a 
prohibitive computational task in large systems. This 
is in stark contrast to kinetic theories based on symme- 
try breaking, where the condensate is a separate quan- 
tity, obeying its own equation of motion. We extend the 
analysis of Ref.|92] to construct an approximate formula 
for the condensate density involving first- and second- 
order correlation functions which are straightforwardly 
obtained. This illustrates the conceptual difference be- 
tween the Penrose-Onsager condensate and the quasi- 
condensate. 

The paper is organized as follows: We first describe in 
Sec. |ll] the procedure of initial state generation within 
each of the two selected stochastic approaches, high- 
lighting the key conceptual differences between them. 
We also briefly review the modified Popov theory that is 
used to benchmark a number of our results. SecHlllad- 
dresses the equilibrium properties such as density pro- 
files and correlation functions, comparing to other, per- 



tinent theoretical results where appropriate. In Sec III D 



the condensate statistics is discussed and some fea- 
tures of the one-body density matrix in the quasi- 
condensate regime are illustrated. Sec.|lV]uses the mod- 
ified Popov approach to address the physical mean- 



ing of the (Penrose-Onsager) condensate mode, which 
may be extracted from the stochastic theories, and con- 
trast this to the quasi-condensate concept. Sec. [V| shows 
that the ncB initial state under GPe evolution does not 
lead to improved predictions for equilibrium properties. 
Sec. |Vl] summarizes our results, with some additional 
material presented in two appendices for completeness. 



II. SIMULATION TECHNIQUES UNDER 
CONSIDERATION 

Each of the stochastic simulation techniques we de- 
scribe here are based on the mapping of a quantum field 
theory of atoms to noisy c-number fields. In this sec- 
tion we discuss them in turn, paying particular atten- 
tion to two important practical elements: (i) the method 
of equilibrium initial state generation and (ii) the nature 
in which the GPe arises as an energy and number con- 
serving means to treat the system dynamics away from 
equilibrium. 



A. Truncated Wigner 

In the truncated Wigner (tW) method, the tempo- 
ral dynamics is governed by the familiar nonlinear 
Schrodinger or Gross-Pitaevskii equation (GPe): 



ih 



dt 



GP[ 



HGp[\m 

2to dz^ 



l/(z)+.9|^(z)p, 



(2) 



where the nonlinear Hamiltonian Hqp contains the 
trapping potential V{z) and the effective two-body in- 
teraction strength g, and /i is the chemical potential. 
While Eq.||2| looks identical to the usual T = GPe 
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equation, the meaning of tjj is quite distinct: (i) The com- 
plex field ip{z,t) represents the condensate, its elemen- 
tary excitations, and the "thermal cloud" surrounding 
it. (ii) The initial conditions are stochastic and include 
quantum and thermal fluctuations of the condensate 
and its excitations. This is essential for incorporating 
spontaneous processes (scattering or decay) that are not 
captured within mean field theory, see, e.g. Ref.| 93ll94l . 
(iii) Averages of operator products are first symmetrized 
before being mapped to classical fields, so that the one- 
body density matrix becomes, for example. 



nq5z 



(3) 



where the average (. . .)vk is taken over the initial condi- 
tions. The second term is a (5-f unction on a spatial grid 
and proportional to the "quantum density" 



(4) 



2Az 



where Az is the grid spacing. As a consequence of 
this "Wigner symmetrization", the density in the non- 
linear term in the Hamiltonian Hq^ [Eq.S] should ap- 



pear with a subtraction, \'4j{z) 



iq. We have incorpo- 



rated the corresponding (small) energy shift guq into the 
chemical potential fi. (See, e.g, Ref.|95| how to general- 
ize the mapping ||3} to two-time correlations.) 

B. Number conserving Bogoliubov initial state 

To obtain a thermal initial state for use in the tW simu- 
lations, we employ a stochastic sampling for the canon- 
ical density operator at thermal equilibrium, based on 
the (number-conserving) Bogoliubov (ncB) approxima- 
tion. In the usual Bogoliubov theory, one shifts the Bose 
field operator by a c-number field (the order parameter), 
which is equivalent to assuming that the system is in 
a superposition of states with different particle number 
(coherent state). The number-conserving version of the 
Bogoliubov theory |[78l[79| . preserves the total number 
of atoms, N, and is constructed to provide the correct 
counting statistics for the condensate mode, in the limit 
of a small thermal component. (The distribution func- 
tion P{Nc) for the number of atoms in the condensate is 
discussed in Sec [III D ) Here, we summarize a practical 
scheme to sample the canonical equilibrium density op- 
erator for the quantum field at a fixed number of atoms 
N, as explained in Refs. |.37, ,46,1 ; a number of technical 
details can be foiind in Appendix [A[ 

1. Condensate mode 

The initial value for the classical field ip{z, 0) is split as 

tp(z,0) = ac0c(z) + 7/'_L(z), (5) 



where the first term describes the condensate, ac being 
the corresponding complex amplitude and Nc = |acp 
the number of condensate atoms. The condensate mode 
function (pd^) is normalized to unity and is given in 
Eq.| |Al| . The splitting ||5j is motivated from an expan- 
sion in the limits of large particle number N and small 
interaction constant g [78 , 80 , 81 , 96] . More precisely, the 
condensate mode and its excitations are calculated self- 
consistently up to second order in (Nth/NY^^, where 
A'th is the number of non-condensed particles, respec- 
tively. (This number also has a small quantum contribu- 
tion.) 



2. Elementary excitations 

The non-condensate field i^_\_{z) in Eq.ljS} arises in the 
next-order contribution of the ncB expansion. It is ex- 
panded in the basis of the Bogoliubov modes: 



k 



[bk Uk{z) + bl vliz)]. 



(6) 



The mode functions are the eigenvectors {uk, Vk)^ to the 
eigenvalue Ek of the Bogoliubov-de Gennes operator 
Cq given in Eq.| |A3| . The Bogoliubov spectrum {Ek} 
is positive and gives the quasi-particle energies relative 
to the chemical potential. The Bogoliubov amplitudes 
bk, b*/. in Eq.||6j are sampled as independent complex 
Gaussian random numbers with zero mean and vari- 
ance (il = (1/2) coth{f3Ek/2) lEZllMl- The mean pop- 
ulation of a Bogoliubov mode is thus equal to 



{\bk\ 



(7) 



the Bose-Einstein occupation number N{Ek) plus an ex- 
tra contribution 1/2. This extra term appears due to the 
symmetric ordering of the quantum operators: 



{\bk 



\')w^l{blbk + b,bl), 



(8) 



and represents quantum fluctuations. Quantum fluc- 
tuations mimic spontaneous scattering into otherwise 
empty modes within a classical field approximation 
L53J . They also lead to the depletion of the condensate 
mode Il97l . The condensate mode function (j)c{z), indeed, 
contains a correction that depends on these amplitudes 
[the field 02(z) discussed after Eq.|A8l]. Within the 
number-conserving scheme of Ref.|38|, this correction 
reflects the change in the condensate number (quantum 
and thermal depletion), the interaction between conden- 
sate and non-condensate particles via the Hartree-Fock 
potential 2g{\4>±{z)\'^)w and via the anomalous average 



.g((i/;_L(^))^)H',seeEqs.|A9 AlOl. 
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3. Condensate number 

The stochastic ensemble of the non-condensate fields 
ip±{z) now determines the sampling of the condensate 
number. The number of condensed atoms Nc is calcu- 
lated from 1^, 38 1 

= N - Nt^{{bk}) + A{{bk}) , (9) 
iVth({6fc}) = AzY,\^±iz)\^ -M/2 . (10) 



Here, A'th gives the number of non-condensed atoms, 
while Ai is the number of terms in the expansion ||6| 
and depends on the computational grid, and acts so as 
to subtract the "one half atom per mode" from the quan- 
tum fluctuations in The quantity A{{bk}), given in 
Eq.( |A7| , averages to zero and implements the "canoni- 
cal constraint" at the level of variances: it ensures that 
the fluctuations of the condensate number Nc are anti- 
correlated to those of A^th (calculated in normal order), 
since the two number operators sum to the fixed to- 
tal particle number N. Once Nc is calculated for each 
member of the ensemble, the condensate amplitude Gc 
in Eq.lS} is taken as Oc — y/Nc- The global phase is ar- 
bitrarily fixed here but the ncB construction is actually 
U(l)-invariant (see endnote [137]), and the phase drops 
out in our observables of interest. A typical example for 
a Wigner field is given in Figjlj 



4. Validity range 

We summarize a few issues that have to be consid- 
ered for the initial state preparation within the ncB ex- 
pansion and the tW scheme. There are two aspects 
here: First, the truncation made in order to obtain the 
evolution equation ||2j assumes that third derivatives of 
the Wigner functional for the quantum field are negligi- 
ble (IZl US |53l. This should be the case when the total 
number of particles is much larger than the number of 
relevant modes in the simulation, i.e. N ^ A4. Given 
this restriction, however, the tW scheme gives approxi- 
mate physical results even beyond the time scale where 
the Bogoliubov theory, in its number-conserving form, 
fails [37J. 

The second aspect is related to the low-temperature 
expansion of the number-conserving Bogoliubov ap- 
proach that is behind the initial state preparation 13711381 
[46 1 . The canonical distribution of the Bose gas is calcu- 
lated by approximating the Hamiltonian of the quantum 
field theory by the quadratic Bogoliubov Hamiltonian 
(as also done in Refs.|71, 98 1). This is valid when the 
number of non-condensate particles is small: Ath ^ N 
which implies relatively low temperatures. The sam- 
pling of the condensate statistics, Eq.||9}, is also an ap- 
proximation that has been discussed in Ref.|37|: in the 
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Figure 1: (Color online) Condensate density Nc\4>c{x)\'^ 
(dashed, black) and a typical realization of the non-condensate 
density \Tj/{x)\'^ = \Tp±{x)\^ (noisy red curve; data is multi- 
plied by 5 to be visible on this scale). The depletion correc- 
tion to the condensate mode, 4>2 (z), is illustrated by plotting 
the 'interference term' 2Re[cl)l{x)(f>2{z)] (dot-dashed, green; 
also multiplied by 5). Note that the condensate atom num- 
ber Nc ~ 19500 depends on the realization. The zero- 
temperature Thomas-Fermi shape for the same parameters (in- 
verted parabola; thin, solid brown line) is also plotted. Here 
(and in most of the paper) we fix fi — 22.41 huj (correspond- 
ing to TV = 20000 for the GPe at T = 0), choosing here 
fcsT — 46 huj. Densities are plotted in imits of g/fi. 



quantum field theory, Nc is restricted to be an integer 
(eigenvalue of the number operator ajflc)/ while the clas- 
sical simulation returns continuous values for Nc- Both 
schemes coincide when the counting statistics P{Nc\ip±) 
(conditioned on a given non-condensate field is 
broad enough, and can be extended to a smooth func- 
tion of Ac. It has been shown that this condition can be 
met when t/j^ is sampled with sufficiently many modes 
{■Mth, say) having a significant thermal occupation (see 
Il38ll for more details). For a typical mode spacing hu, 
this gives a lower limit on the temperature because we 
need A^th = kBT/{huj) ^ 1. Under these conditions, 
one can also justify the gaussian probability distribution 
for the non-condensate field V'i that is used in the sim- 
ulation scheme |38|. 

In summary: since the Bogoliubov approximation is 
used for the representation of the A^-body density op- 
erator, the temperature should be not too high but at 
the same time also not too low because of the smooth- 
ness of the distribution function. This poses limits on 
the applicability of the state preparation protocol within 
tWncB. The above requirements can be checked from the 
inequalities Il37ll 



Xth < (Ath) « <T^{Nc) 



(11) 



where a'^{Nc) is the variance of the (unconditional) con- 
densate statistics P{Nc), obtained from the sampling ||9}. 
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The first inequality precludes very low temperatures be- 
cause A'th would be too small. In addition, the conden- 
sate statistics must not be too broad, 

a(7Ve) « (N,) = iV - (TVth) (12) 

because otherwise the probability of returning negative 
values for Nc would become significant. The method 
gives imphysical results if a large fraction of negative 
values of Nc is returned and we find that this happens 
already at moderate temperatures, when {Nth)/N be- 
comes of the order of ^ 0.2. 

C. Stochastic Gross-Pitaevskii equation 

Within the stochastic Gross-Pitaevskii equation, a fi- 
nite temperature equilibrium state is obtained dynam- 
ically by evolving a complex field tp{z,t) that is cou- 
pled to a thermal cloud which, when approximated as 
in thermal equilibrium, acts as a heat bath (energy and 
particle reservoir) liSllSSllSTllS^ra. 



of a projector operator restricting d3mamics to low en- 
ergy modes, and is termed the stochastic projected GPe 
(spGPe) ESI. 

Technically, the particular sGPe discussed here is ob- 
tained by expanding the system density matrix over 
coherent states using functional integration techniques, 
which leads naturally to the Keldysh non-equilibrium 
formalism. Ultimately, this is found to give a Fokker- 
Planck equation for the Wigner distribution function of 
the entire atomic quantum field II100[[T03 I. This proce- 
dure maps symmetrically ordered correlation functions 
of field operators onto stochastic field correlations, as is 
evident from the fact that each mode k of the classical 
field ip occurs in the stationary limit with an occupation 
number N{Ek) + 1/2 f43!,l60l| (cf. also Eq.0). Within the 
modes which are predominantly classical, i.e. highly oc- 
cupied, N{Ek) ^ 1/2, and this symmetrization is no 
longer important, a common consideration in all classi- 
cal field methods [40-42, 53]. This approximation will 
permit us to move from the quantum relation Eq. |[14]| 
below to its classical counterpart Eq. (TS) which makes 
the simulation scheme much simpler. 



1. System plus bath split 

We may physically motivate a division into two sub- 
systems: the system is represented by the field ijj{z,t) 
and describes the low-lying modes of the ultracold gas. 
These are highly occupied, therefore a classical field de- 
scription is appropriate; the "thermal cloud" of atoms 
whose energies are high above the typical energies of 
the condensate and its excitations, obeys a separate 
quantum Boltzmann equation f43]. Both subsystems 
are naturally coupled to each other by exchanging en- 
ergy and particles, hence the description is given within 
the grand-canonical ensemble. This leads to a nonlin- 
ear Langevin equation (see Eq.||T3|l), often termed the 
stochastic GPe (sGPe). The system dynamics now com- 
bines deterministic aspects (encapsulated within the 
usual GPe) and a stochastic coupling to the heat and par- 
ticle reservoir of thermal atoms. 

We note that there are two distinct formulations of 
such a nonlinear Langevin equation, which are moti- 
vated by the same physical ideas, but which arise from 
very different formalisms (see Ref.[52J for a review and 
a discussion of subtle differences). The derivation of 
Stoof, which we shall adopt in this work, is based on 
the Keldysh non-equilibrium formalism |43 /49"T00'| and 
the resulting theory was first implemented numerically 
in Ref. BSl . An equation that differs in some details 
was formulated by Gardiner, Anglin and Fudge [50 1 and 
cast into its current form by Gardiner and Davis [51J , as 
a limiting case of the quantum kinetic theory put for- 
ward by Gardiner and Zoller 118l ilQllil02 l. The stochas- 
tic equation which results, differs primarily in the use 



2. Stochastic equation of motion 

Making a Hartree-Fock t5^e Ansatz for the probability 
distribution representing the entire trapped gas, leads 
to two separate probability distributions, representing 
separately the high- and low-lying system modes. In- 
tegrating out the low-energy modes, one finds that the 
former may be treated by a quantum Boltzmann equa- 
tion (qBe) [43|. Integrating instead over the high-energy 
modes leads to a nonlinear Langevin equation, 

th^ = (^HgpM^] - - iR{z,t)y + r,{z,t) , (13) 

where R{z, t) is a damping term, which should in gen- 
eral be time dependent, and 77(2, t) a stochastic "force". 
Note that it is for convenience in later discussions only 
that we use the same symbol for the tWncB and sGPe 
wavefunctions, despite the differences in physical con- 
tent. Assuming that the dynamics of the high-energy 
modes may be neglected, the physical picture underly- 
ing this equation is a splitting of the quantum field into 
low-lying modes (the "system", described by ijj{z,t)) 
and a "thermal particle bath", which is considered to 
be at equilibrium and so, on average, Bose-Einstein dis- 
tributed [49,1 (since this is the equilibrium solution to 
the qBE). The term -iR{z,t)ip{z,t) describes the parti- 
cle exchange due to collisions between system and bath 
atoms. The real part of the operator R can be positive 
or negative, corresponding to loss or growth, as for the 
analogous operator appearing in the ZNG scheme [22 J. 
Since collisions occur randomly, the sGPe contains an 
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associated 'noise' term r]{z, t) in Eq.l 13 1. The presence of 
both terms, dissipation and noise, is essential to ensure 
that the fluctuation-dissipation theorem is satisfied: the 
system is thus guaranteed to reach the correct equilib- 
rium at a given temperature. 



3. Damping and noise 



The damping operator is given by the relation 



iR 



hi: 



K 



{N{i,^ + 1/2)- 



(14) 



where KL^ is the so-called Keldysh self-energy (a com- 
plex quantity) and iV(ec — /i) is the Bose-Einstein distri- 
bution, representing the occupation of low-lying modes 
as a function of mode energy e^- While this relation is 
exact at equilibrium, it cannot be easily implemented in 
this form in numerical simulations for the following rea- 
sons: firstly, the mode energy Ic actually corresponds 
to the nonlinear differential operator appearing in the 
GPe |13 i; moreover, is determined by the thermal 
particle distribution, whose accurate temporal represen- 
tation would require solving a qBe. We therefore restrict 
the sGPe to its classical limit, as all current numerical 
applications of this theory do. In the approach of Stoof, 
this means that the action of the damping operator takes 
the simple form 



(15) 



where is still spatially dependent in general ||49l 
l66l , as the thermal particle energies are affected by the 
condensate mean-field. This enables the equation to be 
cast in the form 



ih 



~dt 



(1 - Z7) [HgpM^] + vi^, t) , (16) 



where 7 = ihT.^ /iksT. The form of this equation is 
the same as the stochastic pGPe (spGPe) implemented 
numerically by Davis and collaborators, except for the 
projector Il53l (see Ref . 115211 for a more detailed compari- 
son). 



In Eq.l 16 I, the term ri{z,t) is a complex, Gaussian, 
white-noise process with correlations 



(r/*(z,t)r/(z',t')) = 2h-fkBT5{t - t')5{z - z') 



(17) 



This relation can be read as a fluctuation-dissipation the- 
orem for the system, since the strength of fluctuations is 
proportional to the damping parameter 7, on the one 
hand, and the temperature T of the heat bath, on the 
other This link is essential for the preparation of a ther- 
malized system. A similar stochastic scheme was for- 
mulated in Ref. llSTi for the ideal Bose gas where the 



noise was filtered in order to preserve the total atom 
number (canonical ensemble). 

Since our primary interest in the present work is in 
generating an ensemble of equilibrium states, we make 
a further, numerically convenient, simplification and 
treat as a parameter independent of time and, ad- 
ditionally, space. This is also a standard approximation 
in spGPe simulations [53J, and we have tested that a 
spatially varying does not strongly affect the equi- 
librium state. So, for our present purposes, we solve 
Eq. 1(16) with the dimensionless quantity 7 = 0.01. 



4. Validity range 

The main limitation behind the sGPe used here is 
the classical approximation (highly occupied modes), 
as highlighted by the relation |[15|. In the case of a 
trapped system, this means that the applicability of the 
theory varies spatially. Indeed, the classical approxima- 
tion is better suited to the low energy, central region of 
the trap, in which there are many particles due to the 
presence of a Bose-Einstein condensate. In the outer 
trap regions, there are fewer atoms, and hence there 
comes a point beyond which the classical approach is no 
longer well justified. In general, this point is dependent 
upon the choice of grid spacing, as a finer grid includes 
more high-energy modes. For a given grid with spac- 
ing Az, the accuracy of the classical approximation can 
be checked by comparing, for example, the average den- 
sity (IV'(z)P) to the quantum density' riq = {2lS.z)~^ (this 
value arises from operator symmetrization in Eq.||3} on 
the grid). We shall see that this limits the applicability 
of the classical approximation to the sGPe typically to 
the spatial range where the trapping potential is not too 
large, V{z) — /i < ksT. Within this study, we are inter- 
ested primarily in the central region | z | < i? where the 
condensate is present {R is the Thomas-Fermi radius), 
as also highlighted in the original sGPe numerical im- 
plementation [48] . We have verified that changes in the 
properties which form the basis of our comparison are 
negligible over a range of grid spacings. This is physi- 
cally equivalent to the statement that our equilibrium re- 
sults are unchanged for a range of cutoff energies, which 
mark the split between the low- ('classical') and high- 
lying ('thermal') modes (see also Refs.ll44ll53ll99llT04ll). 



5. State preparation 



We now briefly explain how the sGPe 1 16 1 works in 
practice - more details on this can be found in recent 
reviews Il53l l99l 11041 . As initial condition for the sys- 
tem, one can start with i!{z^ 0) = 0. The dissipative term 



27 in Eq.l 16 leads to a change in the norm of ip{z, t), 
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Figure 2: (Color online) Growth to equilibrium obtained by 
numerical solution to the sGPe. (a) Snapshots of the atomic 
density profile {\ip{z,t)\'^), with time increasing from bottom 
to top. (b) Growth in the average central density vs. time - 
here the times highlighted by colored squares label the corre- 
sponding colored density profiles in (a), (c) Growth in the par- 
ticle number as the system approaches equilibrium. The inset 
shows trajectories from single numerical realizations, which il- 
lustrates the fluctuating particle number between these, com- 
pared to the result of the main plot in (c), which was ob- 
tained by averaging over 1000 such trajectories. (Parameters: 
as stated in subsectionjll F| but with ksT = SGOhcj.) 



but this cannot increase a zero initial condition. It is the 
Langevin force r]{z, t) that 'seeds' the field, as discussed 
in [48]. The particle number N{t) = Az|i/'(z,,)|2 in- 
creases until i}{z, t) relaxes to the solution of the station- 
ary GPe, at a given chemical potential 



Hgp[|V'I']V^ = a^V', 



(18) 



as illustrated in Figj2] The state preparation in the 
sGPe is thus performed dynamically, as the system grows 
to equilibrium in contact with a heat bath at a speci- 
fied temperature. Once the dynamical equilibrium is 
reached, the presence of the noise term ri{z,t) ensures 
that N{t) fluctuates about its final value. The final, aver- 
age atom number (N) depends on the heat bath param- 
eters (temperature T, chemical potential /i), on the trap 
parameters and the atomic species (through the inter- 
action constant g). Although the subsequent dynamical 
evolution requires these noise terms to be maintained, a 
simpler scheme, bearing close analogies to the truncated 
Wigner method, can be based on d5mamical propaga- 
tion of the sGPe equilibrium state via the GPe; such an 
approach was first implemented in 11 641 to discuss quasi- 
condensate growth on an atom chip. 

It is clear that in the grand-canonical ensemble, N has 
a statistical distribution of non-zero width; this width is 
related in the case of the sGPe to the 'history' of the par- 
ticle transfer between system and heat bath as modelled 
by the Langevin seed t). In the simulations, we ob- 
serve indeed that N differs quite substantially from real- 



isation to realisation, however the equilibrium value is 
obtained with reasonable accuracy after averaging over 
a few hundred of them. The temporal variations in the 
ensemble-averaged particle number then become rela- 
tively suppressed. For smoother results and for im- 
proved accuracy, all results presented in this work are 
based on a sample of at least 1000 individual realisa- 
tions. 



D. Linking the theories 

At T = 0, the Gross-Pitaevskii equation represents 
the dynamics of a Bose-Einstein condensate by treat- 
ing the system as made up of a single, coherent mode. 
Within the tW approach, the GPe is instead used to 
propagate a multi-mode system, the initial conditions 
being chosen at random to sample the initial density 
operator of the system. The ensemble of wave func- 
tions ip{x,t) represents all modes of the matter wave 
field, within the limits set numerically by the spatial 
grid, which physically corresponds to the energy cut-off 
choice for the modes being probed. It is for this rea- 
son that no explicit noise terms (Langevin forces) ap- 
pear, and that the total number of atoms (the norm of 
ip{x, t)) is conserved. The initial mode amplitudes com- 
bine quantum and thermal effects [see Eq.|[7||] consistent 
with the Wigner mapping to symmetrically ordered op- 
erator products. Polkovnikov has shown that the trun- 
cated Wigner approximation captures the next order 
correction beyond Gross-Pitaevskii in an expansion in 
h [47] . The initial state that we prepare here is based on 
a fixed number of atoms (canonical ensemble), using 
the number-conserving Bogoliubov approach, although 
alternative (grand-canonical) schemes could be adopted 
as well 1^53 1 . Once we focus on low-lying modes, like the 
condensate mode, we recover nevertheless a broad dis- 
tribution for the non-condensed atoms A'th (the mirror 
image of the condensate statistics P{Nc)). In this per- 
spective, we can even consider /i in Eq.||2j as a "chemi- 
cal potential" for the Bogoliubov modes: the condensate 
plays the role of a particle reservoir, as is perfectly rea- 
sonable if its population is large [ 105 J . 

In the sGPe approach, the wave function ip{z, t) rep- 
resents instead the low-lying modes of the system; al- 
though higher-lying modes should in principle be de- 
scribed by a quantum Boltzmann equation, here they are 
assumed to remain at equilibrium, thereby providing a 
heat bath to the low-energy sub-system under consider- 
ation. In the Bose-condensed phase, low lying system 
modes are highly occupied and the classical approxima- 
tion, amounting to replacing the Bose-Einstein by the 
Rayleigh-Jeans distribution, is well justified. The dy- 
namics of the low-lying modes is quite different, how- 
ever, because particle exchange with the bath is allowed 
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for; this method is therefore a grand-canonical one, as 
can be seen by the growth plots of Fig|2] The "system- 
bath split" can be applied to a trapped gas (a closed 
system) by choosing modes below a suitably small cut- 
off. For a fully self-consistent calculation, in which the 
thermal cloud djmamics are also accounted for, this cut- 
off should be not lower than the global chemical poten- 
tial [43 J. For a classical field method not taking the en- 
tire thermal cloud djmamics into account, this should 
be chosen such that the highest modes simulated are 
macroscopically occupied |53|. It has been proposed 
that a cutoff equal to kgT yields optimum results for 
the condensate statistics of an ideal gas |106|. For the 
purposes of our comparison, we have adopted a differ- 
ent choice here, and have taken for consistency the same 
spatial grid in the sGPe and the tWncB simulations, 
which gives a cutoff of the order of E'max ^ (mAz^). 



Let us summarize the differences between the initial 
state ensembles of the two methods: 

(1) The total atom number N fluctuates in the sGPe 
(grand-canonical), and is fixed in tWncB (canonical). 

(2) The system is thermalized either d5mamically (sGPe) 
by weakly coupling it to a heat bath, or by populating its 
excitation modes with thermal statistics (tWncB). Low- 
lying thermal modes above the condensate equilibrate 
under the sGPe to the Rayleigh-Jeans statistics (classi- 
cal equipartition). This is actually the equilibrium dis- 
tribution for the finite temperature GPe, considered as 
a classical field equation, as has been seen by studying 
thermalization in related classical field methods l.40 -,42J . 
Within the tWncB scheme, these modes are populated 
according to the usual Bose-Einstein statistics, with the 
addition of 1/2 "quantum atom" per mode. This coin- 
cides well with the Rayleigh-Jeans statistics for modes 
with energies below fc^T 1107|, but gives a larger con- 
tribution to high-energy modes, up to the numerical 
cutoff. The tW d5mamics under the GPe redistributes 
these "quantum atoms" with the others, leading, by the 
equipartition law, to an effectively higher temperature. 
This restricts applications of the tW scheme to relatively 
short simulation times. In nearly integrable systems 
(like the quasi-one-dimensional gas), thermalization can 
be quite slow, however II108II109I . 

(3) The energy spectrum of the elementary excita- 
tions is calculated in tWncB approximately, ignoring 
the thermal depletion of the condensate. Indeed, the 
Bogoliubov-de Gennes operator [Eq.|A4l] assumes that 
all particles are in the condensate mode. 



E. Modified Popov theory 

1. Motivation 

While the stochastic approaches discussed thus far 
are suitable to describe both non-equilibrium and static 
properties of the Bose gas, we focus in this study on 
the detailed analysis of a partially condensed Bose gas 
at thermal equilibrium. Mean field theories |52| are 
often applied to study the thermod5mamics in higher- 
dimensional Bose systems, and their solution is, in gen- 
eral, less involved than with stochastic theories. In com- 
paring the equilibrium properties of the sGPe and ncB 
approaches, it will therefore prove useful to have an in- 
dependent method for comparison. 

In lower dimensions, mean field theories have to 
cope with infrared divergences due to the enhanced 
role of fluctuations, as predicted within the Mermin- 
Hohenberg- Wagner theorem UllOillll l. The Popov ap- 
proach f3] where the fluctuations are split into phase 
and density contributions, has proven useful to treat 
phase fluctuations in low dimensions beyond second or- 
der around the mean field. A consistent regularization 
scheme has been developed in the modified Popov the- 
ory of Andersen et al. 171 11121 , extending the work of 
PetroveiflZ. [4] and Kaganef aZ. [5]. The resulting formu- 
las apply to any temperature and dimension, while si- 
multaneously being relatively straightforward to solve. 



2. Quasi-condensate density 

In low dimensions, a condensate does not arise in a 
homogeneous system, but still there is a temperature 
range ^ T < Tc where a so-called quasi-condensate 
can be identified whose density fluctuations are sup- 
pressed lUHZl- Iri the modified Popov theory, the quasi- 
condensate density riqc may be obtained by solving self- 
consistently the following equations for the total density 
[Eq.(4) from Ref.El] 



-y 



N{Ep) 



Er, 



+ - 



^ri Ep 
2Er, 



2ep + 2fi 



(19) 



and for the chemical potential 

H = g{2n - Uqc) 



(20) 



Here, i?p = [ep+2(7?7qcep]^" is the Bogoliubov dispersion 
relation, ep — /2m, V is the system volume, and g de- 
notes the two-body T-matrix (evaluated at — 2/x, which 
corresponds to the energy cost of exciting two atoms 
from the condensate). Eq.|[T9) is evaluated numerically, 
replacing the sums over momenta by an integral. (For 
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Figure 3: (Color online) Illustration of the procedure for solving the modified Popov theory, (left) Graphical determination of 
the self -consistent quasi-condensate density: crossing point between the dashed and solid lines (Ihs and rhs of Eq.|20|, riqc is 
the quasi-condensate density and n the total density). Thin (red) line: classical (high- temperature) approximation to Eq.|T9). If 
we take ^ as in the rest of the paper, the two temperatures correspond to iV ^ 23 800, T ^ 1.3 ^ 0.63 and N = 20 000, 
T « 0.16 ^ 0.074 Tc, respectively (see Eq.jlJ). (right) Temperature-dependent Thomas-Fermi radius -R(r) determined by 
numerically solving for the (local) chemical potential where ii{z) = 2gn[i^{z). The trap is harmonic, chemical potential at the 
center fi = 22.41 hoj as in the rest of the paper, and g = 0.01 (fi^cj/m)^^^. The arrows mark the temperatures chosen for the 
simulations [Table[ll). 



a link with conventional mean field theories see Sec. 
IIIIB2) . 



Eqs. | [19) and pO| | can also be applied in a trap within a 
local density approximation, using a local chemical po- 
tential fi{z) = /i — V{z). The quasi-condensate density is 
then found by solving [Eq.(54) of Ref.[7J] 

(i/Gp["qc] + 2g<h(^)) \/"qc(2) = fi^Juqciz) (21) 

where 2gn[^(z) — 2g{n{z) — nqc(z)) is the Hartree- 
Fock potential due to the non-quasi-condensate parti- 
cles. The spatial point at which 2gn[^^{z) — ii{z) de- 
fines the temperature-dependent Thomas-Fermi radius, 
R{T), that gives an estimate for the thermal depletion 
of the (quasi-)condensate (see Sees. IIIC 1 IV I. We em- 



phasize that this quantity is determined self -consistently 
within the modified Popov theory. For \z\ > B.{T), 
we have nqdz) = 0, and adopting again the local den- 
sity approximation, the atomic density corresponds to a 
thermal gas with a Hartree-Fock interaction. 



n{z) 



dp 



(22) 



e-apip, z) ^ep + V{z) + 2gn{z). 



This procedure is illustrated in FigjSjwhere the left panel 
shows both sides of Eq.|[20|| as a function of the quasi- 



condensate density riqc- The crossing with the dashed 
line determines the self-consistent nqdlJ', T). The (local) 
chemical potential can be lowered, which corresponds 
to moving further from the trap centre, until a minimum 
/^min(r) below which the solution riqc = must be taken 
Il35l . This defines the temperature-dependent Thomas- 
Fermi radius. 



1 - 



i?2 



(23) 



The right panel shows R{T) in the temperature range 
of interest here: the quasi-condensate is shrinking 
smoothly and is about 20% smaller at T ^ T^. 
(This number applies to the parameters introduced in 
SeclHFp 



We recall that the central object within the modified 
Popov theory is the quasi-condensate, but this may be 
linked to the 'true' condensate (if it exists, as defined 
by the Penrose-Onsager criterion), as is discussed in 
Sec. |IV] We will employ this modified Popov scheme 
in order to compare to various equilibrium properties 
where appropriate. This has the advantages of simplic- 
ity and speed over the stochastic approaches, due to the 
relatively straightforward manner in which the above 
equations may be solved. 



F. Parameter choice for comparison 

We wish to address the following issues: 

(i) Initial state generation: how does the finite- 
temperature initial state compare within each method? 

(ii) Ensemble choice: what role does the choice of ther- 
modynamic ensemble play? 

(iii) How does the quasi-condensate and condensate ex- 
tracted from the stochastic approaches compare to anal- 
ogous quantities within the modified Popov theory? 

Our focus being on the relative merits of the stochastic 
methods as thermal field theories, we choose to work in 
a regime in which thermal effects dominate over quan- 
tum effects. We consider a quasi-one-dimensional con- 
finement with a trap frequency oj (oscillator length £) 
and take an effective coupling constant g = 0.01 fkui 
which corresponds to the weakly interacting regime. We 
choose an (average) particle number, iV = 20 000, giv- 
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Figure 4: (Color online) Characteristic temperatures and atom 
numbers of the sGPe and ncB simulations (hollow black cir- 
cles). The filled black circle shows the higher temperature 
regime at which only sGPe simulations were undertaken. The 
(Id) characteristic temperatures in a trap are Tc [Eq.|lJ; Bose- 
Einstein condensation in an ideal gas] and < Tc [Eq.jl); 
phase coherence], shown by the dotted blue and solid black 
lines respectively. The red squares mark the parameters cho- 
sen for the condensate statistics comparison of Sec HID 1| at a 



fixed ratio T/Tc = 0.23, as indicated by the dashed brown line 



ing a chemical potential ^ = 22.41 huj for the ground 
state of the GPe. Within the sGPe, the chemical poten- 
tial is kept fixed as temperature is varied, leading to a 
small variation in the particle number; this is however 
below 6% for the three temperatures at which we un- 
dertake the comparison. These are placed in context in 
Fig|4] In particular, we probe at the lowest temperature a 
regime well suited to tWncB due to the requirement that 
Nth ^ ^/ ari mtermediate regime, and a higher tem- 
perature m which the 'classical' sGPe is expected to be 
most appropriate, due to the occurrence of more highly 
populated modes. These are highlighted in Table|llj also 
showing some other simulation parameters. 





ksT/hu 




T/Tc 


L/R 


M 




(a) low T 


46 


0.052 


0.025 


4.20 


127 


0.22 


(b) intcrmT 


140 


0.16 


0.074 


12.0 


1023 


0.078 


(c) high T 


430 


0.48 


0.23 


17.0 


2047 


0.055 



Table II: Simulation parameters: Atom number and chemical 
potential are fixed to A/" « 20 000, = 22.41 hu, with the char- 
acteristic temperatures Tc and as in Eq.jlJ. The computa- 
tional grid covers z — —L/2 . . . + I//2, with spacing Ajz, num- 
ber of points A^. -R is the Thomas-Fermi radius of Eq.pl) and 
(. the size of the single-particle ground state. 

Length scales in the problem are scaled to the zero 
temperature Thomas-Fermi radius. 



R 



Q.mt 



(24) 



Another relevant length scale is zt = 
\j2{kBT + fi)/fiijj£, which marks the boundary at 
which the thermal energy becomes comparable to the 
trap energy. We take a grid size L > 2zt [Table |ll| 
and a spacing Az/£ = ^27r/(A^ + 1). The number 
of grid points, A4 + 1 = L/Az, is increasing with 
temperature in order to resolve the thermal wavelength 
At = y/27rh/kBT. The grid spacing is chosen such 
that for an ideal gas, the tWncB approach returns the 
correct total particle number (N). In addition, we 
check in the interacting case that doubling Ai does not 
change (N) . For time evolution via the (s)GPe, we use 
a Crank-Nicholson approach, and a fixed time step 
ujAt = 10^^. Ensemble averages are performed over at 
least 1000 noise realizations. 

A typical result of the stochastic methods is given in 
Fig|5] where the data points give the realizations for the 
complex field ijj{z) at selected values of position z and 
temperature. The globally random phase of the sGPe 
data is quite obvious, from the spread around the cir- 
cle, while the ncB ensemble fixes the condensate mode 
to have a real and positive amplitude. The overestuna- 
tion of density fluctuations is also quite visible in the ncB 
data at intermediate and high temperatures, shown by 
the increased variation in the radius of the data points. 
Outside the condensate region (Thomas-Fermi radius 
R), there are no significant differences between either 
method, where the wavefunction represents an incoher- 
ent thermal gas. 



III. EQUILIBRIUM PROPERTIES 

We present here an analysis of the initial states pro- 
duced by the two stochastic methods, that are expected 
to represent thermal equilibrium in the trap. The ac- 
curacy of this equilibrium state is important for mod- 
elling the finite-temperature dynamics when perturba- 
tions take the system far away from equilibrium, e.g., 
changes in the trapping potential. 



A. Density profile 

We begin with the total density profile n{z). Simi- 
lar to experimental data, this contains both the conden- 
sate in the trap center and thermally excited atoms that 
surround it. The equilibrium densities in a harmonic 
trap are plotted in Fig|6] Here and in the following fig- 
ures, the sGPe densities (solid black) are calculated as 
n{z) = {\tp{z,tcq)\'^) where <cq is the preparation time re- 
quired for the system to reach a dynamical equilibrium 
with the bath. In the ncB data (dashed red), we sub- 
tract the "quantum density" rig — l/(2Az) according to 
Eq.| |A12| . This correction is small if the number of atoms 
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Figure 5: (Color online) Data points representing the ensemble of stochastic field values ip{z) at three positions z in the trap 
(from left to right, the average density decreases). Temperature increases from top to bottom. Black: sGPe simulation, red: ncB 
simulation. The dashed green circle indicates, for 1^1 < R, the modulus of the Thomas-Fermi wavef unction, and for z — 1.5 i? 
the square root of the Bose-Einstein density |26} . In the trap center, the Thomas-Fermi approximation yields for these parameters 
|Vtf(0)| « Al/yfL (Atom number N ^ 20 000.) 



per grid cell is large. In addition, its impact on the mean occupation numbers, gives a density 
field is small since gn„ ^ (1 . . . 4) x 10^'' ji. 



For an independent benchmarking, the total density 
profiles are also compared to the total density of the 
modified Popov scheme (green, dot-dashed). At the low 
and intermediate temperature regimes, we find good 
agreement between all three theories. At the higher T, 
Fig|6jc), the ncB result deviates from both the sGPe and 
modified Popov density profiles (which are found to 
agree perfectly with each other [7J). Although the total 
density profiles also include a contribution from ther- 
mal atoms within the condensate region, we defer their 
analysis in this region to Sec III B below. 



Instead, we focus first here on the representation of 
thermally excited atoms outside the condensate region. 
These atoms populate the "wings" i? < Izj < zt, as il- 
lustrated in Figj6|d). In this interval, the gas is still Bose 
degenerate with large occupation numbers per mode. 
Here, the data is well described by a semi-classical ideal 
gas model. At each phase space point {z,p) with ef- 
fective single-particle energy £{p,z) = /2m + V{z), 
assuming the Rayleigh-Jeans law (equipartition) for the 



dp 



2'Kh e{p, z) — pL 



(25) 



V2 T 



arctan 



' E 

V{z) - /X 



Here, E'max = Pmax/2™ is a cutoff energy that depends 
on the maximum kinetic energy on the grid. As shown 
in Figj6|d) (dashed cyan), a good match to the sGPe 
data is obtained for -Emax = 2tt [l / /^.z^hio [ or Pmax — 
2^/tt h/ I^z]. The divergence at z ±i? is an artefact 
due to an infrared divergence of the semiclassical ap- 
proximation. 

Repeating this analysis with the Bose-Einstein distri- 
bution gives a density 



Ap 



2Tih 

ksT 
27r huj 



N{e{p, z) - /i) 



(26) 



.91/2(2;) 



where x = exp[(/i — V{z))/kBT] and the Bose function 
has the asymptotics .91/2(2;) w x for x <^1. Good agree- 
ment with the ncB data is obtained in the limit of infinite 
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5lus quantum 
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Figure 6: (Color online) Average density profiles n{z), normalized as gn{z) / /i, returned by the stochastic Gross-Pitaevskii equa- 
tion (solid black) after an equilibration time tcq and by the number conserving Bogoliubov expansion (dashed red). Plots (a)-(c) 
are for temperatures ksT = 46, 140, 430 /icj listed in Table [ll} R is the Thomas-Fermi radius at T = [Eq.j24|]. The dot-dashed 
green curves give the prediction of the modified Popov theory, as developed by Andersen & al f6l. Plot (d) analyzes, on a loga- 
rithmic scale, the density in the "thermal wings" of plot (c) for R < z < zt where zt (grey vertical line) is the point where the trap 
energy becomes comparable to temperature, V{zt) = ksT + fi. The Rayleigh-Jeans density (dashed cyan) and the Bose-Einstein 
density (dot-dashed green) are calculated for an ideal gas. The brown solid line corresponds to the ncB density plus the quantum 
density level which asymptotes to the latter, Uq [Eq.|4|], indicated by the blue dotted line. 



cutoff, see Fig|6|d) dot-dashed green). Indeed, it is easy 
to see that the contribution of momenta above Pmax is ex- 
ponentially small provided E'max >• fcsT, as is the case 
here. 

At positions beyond zt, the gas enters a non- 
degenerate regime where the occupation numbers are 
small for a large range of momenta. Indeed, the fu- 
gacity satisfies exp[(^ — V{z))/kBT] ^ e^^ for \z\ > 
Zt, and one can make the Boltzmann approximation, 
n{z) cx exp[—V{z)/kBT]. This corresponds to the "ther- 
mal cloud" familiar from the bimodal density distribu- 
tions of a partially condensed Bose gas in a trap |2 1. The 
ncB data provides a smooth crossover into this region 
provided the subtraction of the quantum density is per- 
formed. The brown solid line in Figj6|d) shows the non- 
subtracted density that reduces to a flat background of 
quantum density, n{z) — Uq, in the "Boltzmann tail". 
In this region, the actual number of atoms per mode is 
much smaller than unity, and classical field methods are 
no longer strictly justified. 

By comparing the Rayleigh-Jeans and the Bose- 
Einstein densities [Figj6], it is clear that the sGPe overes- 
timates the number of atoms in the thermal wings. The 
numbers that one gets by mtegrating nR j(z) and ?t.be(^) 
in this region, are quite small, however, when compared 
to the total number of atoms. This is summarized in the 





ks T /hu) 


(iV) [-zt,zt] 


(A'') wing 
RJ-BE 


(iV) tail 
BE 


sGPe 


ncB 


(a) 


46 


20 007 


19 994 


30.9 


7.5 


(b) 


140 


20 132 


19 981 


114.9 


24.6 


(c) 


430 


20 795 


20498 


338.6 


78.5 



last two columns of Table III In the following, we calcu 



Table III: Average atom numbers in the "classical region" out- 
side the condensate and beyond. First and second column: 
average number within \z\ < zt where V{z) < /i+ksT, as ob- 
tained from the numerical simulations (including the Wigner 
subtraction for the tWncB data). Columns three and four are 
based on the ideal gas densities obtained in the classical ap- 
proximation (Rayleigh-Jeans law) and using Bose-Einstein oc- 
cupation numbers. The column "wing" gives the excess atoms 
present in the non-condensate, but still highly populated re- 
gion R < \z\ < Zt outside the condensate (Thomas-Fermi ra- 
dius R). The column "tail" gives the atoms that are located in 
the "tails" zt < \z\ < oo of an ideal gas with the Bose-Einstein 
distribution. For the Bose-Einstein density, an infinite momen- 
tum cutoff is taken, as in Eq.p6). For the Rayleigh-Jeans den- 
sity 1 26 , we take a kinetic energy cutoff -Bmax ~ 2iT{£/Az)^hLo 
with grid spacing Az as given in Table|ll] 



late the total atom number by integrating the density of 
each method between z = —zt-.- zt, where the classi- 
cal approximation is valid. 
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Consider now the significant difference between sGPe 
and ncB in Figj6 c), that occurs within the condensate 
region at the highest temperature. We attribute this to 
the large thermal amplitudes of the Bogoliubov modes 
that are no longer small compared to the condensate 
mode (see FigjS}, and therefore the calculation of the 
non-condensate density has to be done more carefully. 
In the modified Popov theory of Ref.[6J, fluctuations are 
split into contributions to the density and phase, and 
phase fluctuations are systematically discarded when 
calculating the average density. We find that the sGPe 
data agrees well with this total density and are therefore 
confident that it captures the correct result. 



B. Condensate and thermal excitations 

Having considered the total atomic density profiles, 
a further important temperature-dependent quantity is 
the condensate fraction. The classical field methods un- 
der consideration here provide a unified description of 
the lowest modes of a trapped Bose gas, and so further 
analysis is necessary to isolate the condensate fraction, 
as in experiment. We choose to compare the coherent 
and incoherent phases of the gas, and so focus on a 
partitioning based on those atoms within the Penrose- 
Onsager ground state, and those in orthogonal states 



1. Density profiles and depletion 

The phase coherent fraction of the gas is identified 
by making use of the Penrose-Onsager (PO) criterion 
for Bose-Einstern condensation |113|. The stochastic 
simulations provide us, again through the operator- 
classical field correspondence, with the so-called one- 
body-density matrix z') — {ijj* {z)'tp{z')) . This ma- 
trix is hermitian and positive; evaluating its trace by 
integrating spatially, we get {N) in the sGPe method, 
and N + M./2 within the ncB approach. The PO crite- 
rion states that Bose-Einstein condensation has occurred 
when the largest eigenvalue of p{z, z'), denoted here by 
{Nc), is comparable to N , the other eigenvalues being 
much smaller The eigenvector corresponding to 
{Nc) gives us the (PO) 'condensate mode', whose spatial 
width characterizes the long-range phase coherence of 
the degenerate Bose gas. This mode provides a numer- 
ical way to implement the splitting in Eq.ljS} between 
condensate and excitations, since we get the condensate 
amplitude from 



Figure [tI plots the condensate density nc{z) 
{Nc)\(j)c{z)Y' and the 'thermal density' nth{z) ~ n{z) 
ndz), with corresponding thermal fractions given in 
Table 



IV 



ttc = Az^0*(z)?/;(2 



(27) 



and ^j_{z) := ■0(z) — ac4>c{z) is by construction orthogo- 
nal to 6c- 



At the lowest temperature (left images), the 
condensate dominates [Figl7|a)] and the thermal den- 
sity is globally weak [FiglTTd)]. The two broad peaks 
at z ^ ±i? arise because the repulsive interaction with 
the condensate pushes the thermal component out of the 
trap center. In the intermediate temperature data (mid- 
dle images), the thermal fraction is larger, yet the two 
methods give good agreement, with only a marginal 
difference in the peak value of the condensate density 
which now becomes noticeably smaller than the T = 
solution for the same /i (thin blue line). At the rela- 
tively higher temperature fcsT — A30huj (right images), 
the two methods disagree: their condensates are sim- 
ilar in axial extent, but the tWncB mode (dashed red) 
has a lower peak than the sGPe(solid black), so con- 
tains fewer atoms, and as a result the thermal fraction is 
higher. For all temperatures, the approximate conden- 
sate mode constructed in the ncB theory [4>c{z) as given 
in Eq.ljATJ; dot-dashed, blue] agrees well with the PO 
condensate extracted from the density matrix of the ncB 
simulations (dashed red): we show data in Fig|7|[c) only 
as these quantities become indistinguishable at lower 
temperatures. 

Although the modified Popov theory inherently 
solves for the quasi-condensate density, nqc(z) (rather 
than the phase coherent, Penrose-Onsager condensate 
mode plotted here), we also show in Fig. [7 c) a pre- 
diction for the phase coherent condensate which may 
be extracted from the modified Popov approach (dot- 
dashed, green curves), and whose calculation does not 
require the full one-body density matrix, as we discuss 
in more detail in Sec |IV| At this moderate temperature, 
this prediction agrees well with the sGPe, except for a 
small region near the centre; at lower temperatures, we 
have found an even better agreement. 

The difference between the ncB and sGPe results is 
likely due to the overestimation of the thermal density 
in the condensate region within ncB. Consider the com- 
plex values of i^{z) in the bottom images of FigjS] for 
z < R: the density is determined by the modulus of 
tACz), and the tWncB data clearly have more points with 
larger modulus. Another reason may be the way the 
Bogoliubov energy spectrum is calculated: indeed, it 
is based on a condensate wave function, denoted (po in 
Eq.| |Al| , which contains all the particles of the system, 
whereas an improved approach would account for the 
depletion of the condensate in calculating the Bogoli- 
ubov spectrum. To illustrate this, we consider as a first 
step an improved spectrum based upon a a condensate 
with the same number of atoms as contained within the 
sGPe PO condensate. The thermal density which results 
is shown as the dotted maroon line in Fig. [7|f ), and is al- 
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condensate density 
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Figure 7: (Color online) Condensate (top row) and thermal cloud (bottom row) densities for temperatures (a),(d) fcsT — 46fta;; 
(b),(e) fcsT — UOhco; and (c),(f) ksT = 430/ia;. Condensate density nc{z) = {Nc)\cj)c{z)\^ obtained from Penrose-Onsager 
analysis of the one-body density matrix [sGPe: solid, black; tWncB: dashed, red] . Also shown for reference is the T = stationary 
solution to the GP equation with eigenvalue /i (thin solid, turquoise), which coincides with the zeroth order condensate mode 
within ncB. In (c) the condensate density (Nc) 1 02 (2) P following the 2nd order correction within the ncB expansion is also shown 
(dot-dashed, blue), as well as the modified Popov condensate density (dot-double dashed, green) [see Sec|IV|Eq.l|56} for details]. 
Bottom row: thermal density nth(z) = n{z) — ndz) [sCPe: solid black; tWncB: dashed red]. As in Fig|6 the Wigner correction 
was made for the tWncB case. Shown in (d-f) is the T = Bogoliubov prediction, Eq.|A4|, with A'^ (solid brown) and, in (f) only, 
also A'^ — Npo of the sGPe (dotted maroon), alongside the modified Popov result (dot-double dashed, green). 





ksT /hu] 


sGPe 


tWncB 


(a) 


46 


0.0474 


0.0468 


(b) 


140 


0.141 


0.143 


(c) 


430 


0.365 


0.450 



Table IV: Thermal fraction (A'^th) / {N) versus temperature for 
the sGPe and tWncB initial states for the three chosen temper- 
atures, where (A''th) is the integral over nth(z) in the region 
2 1 < zt, and the correction due to symmetric operator order- 
ing is applied to the tWncB data as for Fig[6] 



ready in better agreement with both the sGPe and modi- 
fied Popov results, despite the quite rudimentary nature 
of the modification to the Bogoliubov spectrum. 

The thermal fraction, (iVth)/(^> = 1 - {Nc)/{N), 
varies slightly depending on the simulation method, as 
shown in TablelTVl Since the total atom number N varies 
across the statistical ensembles, we calculate the conden- 
sate fraction from {Nc)/{N). We observe again larger 
values for this quantity in the tWncB method at higher 
temperatures (row (c)). We suggest that the depletion in 
this range is not small enough to warrant the ncB expan- 
sion aroimd a "large" condensate. Indeed, if we identify 
from Eq.(14), Ref.|80j, {Nth/Nc)^^^ as a small expansion 
parameter, within the ncB calculation this reaches the 



value « 0.90 that is clearly not small. 

2. Condensate shape 

The back action of the non-condensate particles can 
be made visible by a careful analysis of the shape of the 
condensate wave function 4>c{z), the results of which 
are summarized in FigjS] The simplest generalization 
of the GPe that applies to nonzero temperature is a 
Hartree-Fock (HP) potential due to non-condensate par- 
ticles (analogous to Eq.|[2T)) 

HP: (iJGp[A^c|0cP]+2.9nth(2))0c = A^0cW (28) 
where the thermal density is 

n,^{z) = (V31(z)7^i(z)) (29) 

One could however also take into account the anoma- 
lous average m(z) due to non-condensate modes. If the 
condensate field ac.(f)c is chosen real, the anomalous av- 
erage is simply given by 



(30) 



(A definition not based on U(l) symmetry breaking can 
be found in Eq.|[47| below.) Within the Bogoliubov ap- 
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Figure 8: (Color online) Analysis of the back-action of the non- 
condensate atoms on the shape of the condensate mode ^c- 
We plot the effective potential of Eq.pi} with (lower curves, 
HFB theory) and without the last term (upper curves, HF the- 
ory), normalized to fi. The condensate mode function <j)ciz) is 
taken from a Penrose-Onsager analysis of the one-body den- 
sity matrix. Solid black/brown lines: sGPe data, (dot-)dashed 
red/green lines: tWncB data. The HFB theory that includes 
the back action from the anomalous average is closer to the 
actual chemical potential. 



proximation, one has 

nth{z) - 5]{iV(i?fe)|ufe(z)|2 + (7V(i?fe) + l)|i;fe(z)|2}, (31) 
fe 

m{z) = Y.(^N{Ek) + l)uk(z)vl{z). (32) 
fe 

(In a homogeneous gas of dimensions D > 2, the sums 
in Eqs.(31 32 1 are ultraviolet divergent ITOl ITTl l85l and 
are regularized routinely by a renormalized coupling 
constant. This is not needed in the one-dimensional case 
considered here, due to the decay of Vk{x) for S> ^.) 
This leads to the following 'Hartree-Fock-Bogoliubov 
(HFB) extension of Eq.||28): 



HFB: 



{HGp[N,\(l>,\^] + 2gnthiz)) (f), 
+ gm{z)(j)l 



(33) 



where we allowed momentarily for a complex-valued 
condensate wave function to illustrate the U{1) invari- 
ance of the theory |137J. 



The data shown in Fig|8] taken from both stochas- 
tic methods, suggests that Eq.(j33| is more appropriate 
for the (PO) condensate mode, at least in the central 
region of the trap. This effect, along with the role of 
higher anomalous averages, has been studied in de- 
tail in the context of the microcanonical pGPe theory in 
||Tl4J. There are significant differences at higher temper- 
atures in the tWncB data, similar to those appearing in 
the average density. We recall that the HFB theory has 
been put into question because it leads, for a homoge- 
neous system, to a gapped excitation spectrum, in con- 
tradiction to the Hugenholtz-Pines theorem |30-34, 52|. 

One should also compare to the modified Popov the- 
ory (Sec II E I. Indeed, we can interpret the HFB potential 
in Eq.||33[l for the condensate mode, as a modification of 
the thermal density. Going back to a real-valued con- 
densate 4)c{z), the effective potential in Eq.|33 1 takes the 
form 

VnFB = V + gN^cfl + 2gntY, + gm (34) 

while the last term is missing in Eq.|[28). Now in Bo- 
goliubov theory, we have (adopting real mode functions 
normalized according to Eq.||A6)) 



+ {2Nk + l)uk{z)vk{z)} 
^ ^{Nk{uk{z) + Vk{z)f + {uk{z) + Vk{z))vk{z)} (35) 

k 

with the shorthand Nk = N{Ek). Equation |35) is 
closely related to the thermal density n^j^ of the modified 
Popov theory [Eq.|[l9|] which can be expressed in terms 
of Bogoliubov amplitudes (normalized to — — 1) 
as 

1 



'th 



= n — = 



V ^ 

p 



{Np{up + Vp) 



+ {Up + Vp)Vp + 



2ep + 2/i 



(36) 



} 



This clearly contains an additional term; this term is, 
however, small in the temperature regime considered 
for our numerical simulations, leading to a good anal- 
ogy when the above re-interpretation of mean field po- 
tentials is made. This is in fact somewhat analogous to 
the work of Ref. lTTS] . 

The analogy cannot be pushed further, since the mod- 
ified Popov theory does not directly deal with the con- 
densate mode (in the PO sense), but with the quasi- 
condensate. See Sec IV for a link between the two quan- 
tities. 



C. Correlation functions 

The focus in this section is upon spatial coherence in 
phase and density. These show a rich physics in weakly 
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Figure 9: (Color online) First order correlation function g*^' (z) — g'^' (0, z), as defined in Eqj 37 <, from the sGPe data (solid black) 
and the ncB data (dashed red). Temperatures in (a-c) as given in Table in (d), we have T ^ 1.3 Tj, (only sGPe data shown). We 
plot the real part of g*^^' (2), the imaginary part being negligibly small. Dot-dashed green lines: Eq.i 39 , based on modified Popov 
theory. Brown, dotted lines: Eq.|40), In (d), we plot Eq.|40) using R{T) < R (solid thin brown) and R{T) — R (dotted orange). 
The contributions from the sGPe Penrose-Onsager condensate mode (dashed blue) and from non-condensate modes (dot-dashed 
red) are shown separately, from Eq.(43). In (a-d) the vertical dashed thin lines indicate R{T). (e) Estimated temperature based on 
fitting g^^'' (z) near 2 = to Eq.(42| (dashed line) for sGPe (black circles) and ncB (red squares). 



interacting one-dimensional Bose systems due to the 
separate characteristic temperatures for the suppression 
of phase and density fluctuations. 



1. First-order coherence: phase fluctuations. 

To study the phase coherence, we begin with the first- 
order coherence function (a normalized one-body den- 
sity matrix) 



the Ansatz 



[niz)n{z')]^/^ 



(37) 



where n{z) is the average density and the normalization 
gives g^' ^^iz, z) — 1. For simplicity, we fix one position 
in the trap centre, z' = 0, and write g^^\z) = ^'^•'(O; z). 
As illustrated in Fig|9| g^^^ [z) scales roughly linearly in 
the centre and drops quickly to zero towards the border 
of the condensate mode (z ^ R). The slope in the cen- 
tre agrees well between the two methods (sGPe vs. ncB: 
solid black vs. dashed red), but for z ^ i?/2, differences 
appear at the highest temperature. Quite striking are 
the negative values for g^^' (z) which then occur, within 
the ncB formalism, in the region of the Thomas-Fermi 
radius (anti-correlation between central region and con- 
densate edge). 

At the low temperatures probed, the behaviour of 
^'^^(z) compares well with the theory of phase fluctu- 
ations in weakly interacting Bose gases: One starts from 



gW(z)=exp[-((^(z)-^(0)f)/2] 



(38) 



and works out the thermal statistics for the phase oper- 
ator 6{z). 

From this point, the required exponent may be cal- 
culated within the modified Popov theory ll7l lll2IITT6l . 
Making the classical approximation N{Ej) « ksT/Ej, 
we can write this exponent as 



{W)-em^) 

AT 



2j 



E 

i>o 



3{j + 1) 



^ (P,{z/R{T))~P,{Q)f 



P,{z/R{T)) 
l~z^/R(Ty 



(39) 



where Pj{z) are Legendre polynomials of order j. Here 
ksT^ « 40 /i is the characteristic temperature for phase 
fluctuations [Eq.|[TJ]. 

The first term of Eq.|p9| was derived by Petrov et 
al. |4] who, focussing on the temperature range hu ^ 
fc^T, did not explicitly include in their expression the 
temperature dependence of the Thomas-Fermi radius 
R{T). For the parameters considered, the second term 
within the sum of Eq.|p9) typically gives a small contri- 
bution, but leads to a rounding off of the central peak 
in ^'^^(O, z). This additional term leads to a divergence 
in g(^) (0, z), however, as \z\ R{T), due to the assump- 
tion that the condensate density is parabolic. The results 
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of Eq.|[39) are shown by the dot-dashed, green curves in 
Fig. |9] In particular, at the highest temperature (panel 
(d)), the shrinking of the quasi-condensate {R{T) < R) 
is clearly seen, due to the increased thermal component 
of the density. 

A closed form for g^^\0,z) may be obtained by ne- 
glecting the second term of Eq.|[39|, since the mode sum- 
mation for the first term can be performed analytically, 
as pointed out previously in [4, 117J. This gives 



mz)-e{o)? 



4T 



log 



l + z/R{T) \ 
\-z/R{T)) 



(40) 



where we have kept the temperature-dependent 
Thomas-Fermi radius, thereby generalizing the expres- 
sions of Refs. ill [1171 . We note that Ghosh considered a 
quantized hydrod5mamic approach to the correlations 
of a quasi-ld Bose gas II118I , which also extends the 
work of Petrov et al. in Ref. flj. Eq.|[40| is shown with 
R{T) taken from the modified Popov simulations by 
the solid, thin brown lines in Figs. |9|a-d) and agrees 
well with the numerical data. At the highest temper- 
ature [Fig. |9];d)], we see a clear deviation between an 
approach like the modified Popov theory that takes 
condensate depletion into account (taking R{T) < R, 
solid, thin brown), and the original expression of 
Ref. [4 J where the condensate size is fixed at a constant 
Thomas-Fermi radius (dotted orange). 

In the central trap region \z\ <^ R, we can approxi- 
mate Eq.|[40) as 



mz)-9m' 



2z 



L^Ty 



AT 



R{T), (41) 



where L^{T) is the phase correlation length. This ex- 
pression illustrates that as T > T^, the system is (first- 
order) coherent over a scale significantly shorter than 
the Thomas-Fermi radius lH. Assuming R{T) « R 
and using the Thomas-Fermi formula for the conden- 
sate profile, we can re- write the product T^i? in Eq.<|4T| 
in a model-independent way 



1^1 g'-'Hz) 



9 



z k^T 



2fiuj£ 



(42) 



Based upon Eq.||42]|, the slope of g^^'^iz) in the region 
< z < R/2 yields an independent measure of the 
temperature. We compare the temperature extracted in 
this way (denoted Tfn) to the input temperature of the 
simulations. Both simulations give coherence functions 



g^^^(2:) that are consistent with Eq.|42 i, except for the 



ncB data at T sa 0.25 where the phase coherence is de- 
caying faster. The ensemble prepared by ncB in the latter 
case appears not only to be at a higher temperature (as 
in Figj7 c)), it is actually not even stationary, as we illus- 
trate in Sec|Vl Let us come back to the link between the 
phase coherence function g'^^ and the one-body density 



matrix. As explained in Sec III B} the latter can be ex- 
panded in orthogonal modes, with the (PO) condensate 
mode (l>c{z) giving the dominant contribution. In the no- 
tation of the ncB approach [Eq.||5|], one can decompose 
the stochastic wavefunction in condensate and excita- 
tion parts, ac(j)c{z) + ip±{z), where ip±{z) represents all 
modes orthogonal to the PO mode. The first order cor- 
relation function may then be broken down into 

irmiz)) ={N,)^:{o)Mz) + (v4(o)v^i(^)>. (43) 

This gives two contributions to g^^^z) that are illus- 
trated in Figj9|d). The condensate mode (dashed blue) 
provides the largest contribution and is positively corre- 
lated. The contribution due to modes orthogonal to the 
PO mode (dot-dashed red) becomes negative towards 
the condensate border because these modes have ad- 
ditional nodes (only even mode functions contribute). 
This reduces g'^'(z) below the condensate contribution. 
The "spike" near the centre is also due to higher modes 
and contains the approximately exponential decay due 
to phase fluctuations [Eqs.|38 |4Tl] that becomes nar- 
rower CLST Ta,. 



2. Second-order coherence: density fluctuations 

We now consider correlations of order four, namely 
fluctuations of the atomic density. These are captured 
by the coherence function 



[n{zW 



(44) 



In terms of field operators, we actually consider the 
probability of detecting two atoms at z, g^'^^{z) cx 
(J\j/t^t\j;\ii(^ This operator ordering is mapped to the fol- 
lowing combination of tW data: 

(^-t^t^^) ^ _ 4(|V'(z)P)w", + 2nl (45) 

where Uq is the quantum density level on the grid 
[Eq.||4|]. The local density n{z) used for normalization 
is also Wigner-corrected and given by Eq. | |A12 1. The 
sGPe data are taken as in Eq. | [44| . 

It is well known that for a single-mode coherent field, 
^(^^(z) = 1, while for a chaotic (multi-mode) field 
with Gaussian statistics, g^'^\z) = 2 IITT91IT20J . Anti- 



bunching, g'^^ (z) < 1, is a non-classical effect that we 
do not expect within the classical stochastic theories 
used here; it occurs indeed at lower temperatures, see 
Refs. lTSl 11211 . Any value in between the limits 1 and 2 
is thus a measure of how many modes effectively con- 
tribute to the density. The data shown in Fig|TO| follows 
the expected behaviour Il44l [1161 [1201 IT22ll : a flat plateau 
in the center and a jump from 1 to 2 at the border of 
the condensate. (The oscillations outside the center are 
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Figure 10: (Color online) Density correlation function ,g^^' (2) [Eg .' 44 ] from the sGPe simulations (solid black) and the ncB simu- 
lations (dashed red). For the ncB case, the corrections of Eq.i 45 are applied, so that in the quantum field theory, g'^' (z) has the 
meaning of a second-order coherence function. The vertical dashed thin lines indicate R{T) at the temperatures in (a-c), which 
are as in Table|llj (d) Comparison of results for g'-'^^ (0) (black circles: sGPe, red squares: tWncB) with Lieb-Liniger theory, Eq.i 46 , 
taken from Ref.|75| (dashed brown). 



statistical errors that are enhanced by the normalization 
in Eq.|[44|.) The jump at the Thomas-Fermi radius be- 
comes more gradual as the temperature rises, and the 
single-mode region shrinks [116J. At the highest tem- 
perature, the ncB theory gives anomalously large values 
g^^-* (z) > 2 near the condensate border. 

Since the lowest excitation modes of the condensate 
carry mainly phase fluctuations 12|, we expect a signifi- 
cant deviation from g^^^ (z) = 1 to set in at a higher tem- 
perature compared to ^(^'(z). This can be made more 
precise by comparing to Ref. llT^ where Kheruntsyan et 
al. use exact solutions of the Lieb-Liniger model, within 
the local-density approximation, to calculate the density 
correlation in a trapped gas. Their result for the trap cen- 
ter, within the weakly interacting regime /i <C 2k bT <C 
keTd = Nhujoi interest here, is (Eq. (5.10) of L75J ): 



g(2)(0)«l + 



4\/2T 
3Td ■ 



(46) 



The linear increase in temperature is in good agreement 
with the results of our classical field simulations, see 



Figure 10 d), and fTT6l. The high-temperature ncB data, 
however, is much too large compared to both the sGPe 
and Lieb-Liniger theory. This suggests either a higher 
temperature, consistent with the findings of previous 
tests, or a significant overestimation of density fluctu- 
ations. We mention that values g^^^ {z) — > 3 II123I would 
arise when the field ijj{z) has a fixed phase and be- 
haves like a real-valued random number. This is related 
to a large contribution from the 'squeezing correlation' 
m{z) ~ {%Ij{z)%1^{z)) that we discuss now. 



3. Squeezing and anomalous average 

We finally consider the anomalous average of the non- 
condensate field defined as (cf. fl24]) 



m{z) 



[{ac4'c[z))*'ip±{z)] 
\ac(t)c{z)\'^ 



(47) 



where ac4>c{z) is the component of the matter wave field 
along the condensate mode [Eq.|[27|] and ip_\_{z) is the 
perpendicular component. Note that m(z) as defined 
in Eq.||47| is invariant under global phase transforma- 
tions of both the total field ipiz) and the condensate 
mode function (j)c{z). It vanishes if the condensate am- 
plitude Oc and V'-l(^) have no fixed phase relation: we 
thus probe the phase locking between the condensate 
and non-condensate fields. We use for our data anal- 
ysis the Penrose-Onsager condensate mode introduced 
in Sec III B This interpretation of the anomalous average 
can be re-phrased in the squeezing language of quan- 
tum optics L125J: the interference term between conden- 
sate and non-condensate fields is split in two quadrature 
fields (both have the dimension of a density) 



{ac<t)c{z))*''J^±{z) Xniz) + iXg{z). 



(48) 



The real part X„ indeed gives the (local) density fluctu- 
ation on top of the condensate density |aci?!>c(^)P, while 
the imaginary part Xg describes phase fluctuations (if 
these are small). On average, these quadratures are zero, 
and the difference of their variances is 



{X^{z)) - {Xliz)) = Re([(a,0,(z))>i(z) 



(49) 



which is just the real part of m(z) defined in Eq.|47 1. The 
sum of these variances equals the normalization factor 
in Eq.|47| times the non-condensate density nt\i{z) = 

Within the Bogoliubov approximation, we calculate 
m{z) by choosing a phase reference where the conden- 
sate field (pdz) is real-valued. By expanding '^l)^_{z) 
over Bogoliubov modes with operator amplitudes (5k — 
a*hk/\ac\ instead of hk, global phase invariance holds 
(see Refs.l70l l80l[T37 l). In the Bogoliubov limit, conden- 
sate number fluctuations can be ignored, jocp — {Nc), 
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Figure 11: (Color online) Real part of the squeezing correlation (anomalous average) m{z), as defined in Eqjj47|, for the tempera- 
tures of Table Solid black: sGPe result, dashed red: ncB result. The Bogoliubov result (solid brown, Eq.|32^) is calculated with 
the mode functions of the T = BdG operator jA3} , and also using the Bogoliubov spectrum calculated for a condensate number 
equal to the sGPe PO condensate number (blue, dot-dashed). 



and we recover Eq.|[32) 

m{z) = ^(K(z)/3,+i;;^(z)/3^)K(z)/3,-|-«^(z)/?i.)) 
fe 

= Y,{'^N{Ek) + l)uk{z)vl{z) (50) 

k 

where TV = {(31(3^,). This quantity is thus sensitive 
to the 'anomalous' or 'hole' part Vk{z) of the Bogoliubov 
modes. In particular, we note that the anomalous av- 
erage shows a quite precise linear scaling in T in the 
temperature range of interest. This illustrates the rela- 
tive dominance of highly populated modes that are well 
described within the classical approximation. 

The data in Fig 11 show a reasonable qualitative 
agreement between the stochastic simulations and Bo- 
goliubov theory: m(z) has a large negative real part. 
Beyond this, there are clearly differences on a quantita- 
tive level, particularly as temperature is increased. Note 
that the anomalous average is of comparable magnitude 
to the non-condensate density nth{z), which points to- 
wards a strong enhancement of the phase fluctuation 
quadrature Xg{z) relative to density fluctuations. The 
agreement between sGPe (solid black) and ncB (dashed 
red) remains reasonable at all temperatures considered 
in our comparison, however both theories deviate from 
the (T = 0) Bogoliubov prediction (brown, solid line). 
Also shown in Fig. [ll]is the result of Eq.|[50), calculated 
with Bogoliubov mode functions for a condensate with 
a reduced number of atoms: we have replaced in the 
Bogoliubov-de Gennes operator iV|0oP by {Nc)\(f>c\'^, 
where (Nc) is obtained from the sGPe. Note that this 
does not significantly improve the agreement, unlike the 
case of the thermal density of Figl7|f). The enhancement 



of phase fluctuations is nicely illustrated in Fig 12 where 
the realizations of the complex field tp{z)a*/\ac\ are plot- 
ted. The additional phase factor, relative to the data 
of Fig|5| removes the random phase of the condensate 
mode, and reveals the squeezed distribution of the com- 
plex field, with enhanced fluctuations in the quadrature 



orthogonal to the condensate mode. At higher temper- 
atures, the sGPe data show that these fluctuations are 
channeled into a "crescenf'-shaped region, maintaining 
the relative suppression of density fluctuations. The ncB 
expansion does not take this into account, and the fluc- 
tuations keep their alignment to orthogonal quadratures 
so that density fluctuations (in the radial direction) be- 
come too large. This clearly happens when the phase 
difference across different points in the system becomes 
comparable (in standard deviation) to tt/2 so that lin- 
earization procedures break down ("quasi-condensate 
regime"). 



D. Condensate statistics and fragmentation 

We analyze in this section the one-body density ma- 
trix (z)V'(z') in more detail. Its eigenvector with the 
largest eigenvalue corresponds to the condensate mode 
(pcix) in the sense of Penrose and Onsager, as explained 
in Sec IIIB| The distribution function of the correspond- 
ing complex amplitude ac [Eq.pTjl] provides us with 
the probability of finding Nc ~ |acp atoms in the con- 
densate, the so-called "counting statistics". We empha- 
size that this quantity depends on moments (correla- 
tion functions) of arbitrarily high order of the stochastic 
field. We also discuss the relative importance of non- 
condensate modes whose occupation grows as the tem- 
perature increases, illustrating a phenomenon similar to 
fragmentation HI for T ^ and above. 



1. Counting statistics 

The statistics PiN^) of the number of condensate 
atoms has been well studied in the context of laser the- 
ory [24J, and Bose-Einstein condensation ll7TU73ll83ll84l 
l98l 11261 1127II . It is worth mentioning that number dis- 
tributions for an ideal Bose gas provide an example 
where the canonical and grand-canonical ensembles of 
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Figure 12: (Color online) Representation of the ensemble of stochastic fields ipiz), locked to the phase of the condensate mode (we 
multiply with a*/\ac\ where Oc is the amplitude l |27| of the (Penrose-Onsager) condensate mode). Black: sGPe simulation, red: 
ncB simulation. Trap position and temperatures as in Fig|5] Note the elliptical shape of the point cloud, illustrating the relative 
enhancement of phase fluctuations (imaginary part). The sGPe data deform into a "crescent" shaped cloud at intermediate and 
high temperatures, while the ncB distribution remains aligned to orthogonal quadratures. The dashed green circles represent the 
Thomas-Fermi {\z\ < R) and Bose-Einstein {z — 1.5 R) predictions for the square root of the density at each trap position and 
temperature, as in Fig|5] 



thermodynamics are not equivalent: only the average 
atom numbers coincide, while all higher moments of 
the (total) atom number are anomalously large in the 
grand-canonical ensemble 1 70 , 84 1 . This anomaly is re- 
moved in an interacting gas due to the energetic cost of 
adding particles to the condensate. The counting statis- 
tics P{Nc) is found from the stochastic data by drafting 
a histogram of the values for Nc — |acp across the en- 
semble of realizations, with Oc calculated from Eq.|[27|. 
Obviously, the Nc need not be integers here, due to 
the replacement of operators with classical fields. Fig- 
ure 13 compares these data with the theory of Scully 



and co-workers (S&Co), in particular Ref.|71J. This is 
developed within the canonical ensemble, and treats the 
non-condensate modes either within Bogoliubov theory 
(low temperatures), or extrapolated to higher tempera- 
ture with the help of a rate equation approach. We out- 
line the main steps in Appendix |B] 

Reasonable agreement between the stochastic simula- 
tion methods and S&Co (dot-dashed green) is apparent 
at all three temperatures shown in Fig 13 At the lowest 
temperature [Fig 13 a)], the sGPe approach (solid black) 
gives a slightly a broader distribution. This broaden- 



ing becomes more pronounced if the atom number is 
lowered [panel (d)], bringing the system closer to an 
ideal gas. At the highest temperature considered, the 
ncB distribution (dashed red) gives too much weight on 
small condensate numbers. This failure is particularly 
striking, given that the method of S&Co uses a Bogoli- 
ubov description of the non-condensate particles that 
is fairly close to the ncB expansion. There is one addi- 
tional ingredient, however, namely the growth and de- 
pletion rates in the rate equation Ansatz for the counting 
statistics (see App. B for details): these rates are calcu- 
lated as a function of Nc, while in the ncB expansion, 
the Bogoliubov spectrum is calculated only for the ex- 
treme case Nc = A^. We thus expect that some effects of 
a strongly depleted condensate are not captured. This 
illustrates again the importance of self-consistently ad- 
justing Nc within the theory as temperature is varied, 
see also Ref.ll72l. 
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Figure 13: (Color online) Counting statistics P{Nc) for the condensate number, with Nc — l^cp, and the condensate mode in 
Eq.|27| obtained from the Penrose-Onsager criterion. Solid black: sGPe, dashed red: ncB, dot-dashed green: theory of Ref.|71|. 
Temperatures in (a-c) increase from left to right, (as in Table|llf, (N ^ 20 000, ^ = 22 Al hoj). Panel (d) has TV ^ 1 000 (/.i = 3.11 huj) 
and a temperature T « 0.23 Tc, the same ratio as in (c). Dotted blue line: ncB with Rayleigh-Jeans instead of Bose-Einstein 
occupation numbers, i.e., in Eq.|7|, al ^ keT/Ek. 



2. Discussion 

It may come as a surprise that a grand-canonical ap- 
proach like the sGPe where the total atom number is not 
fixed (i.e., values Nc > (N) are not excluded), is able to 
reproduce the counting statistics of number-conserving 
theories (like the ncB and S&Co). We attribute this to 
the interatomic interactions in the system that trans- 
late fluctuations in the condensate number into ener- 
getic changes. This makes the system "stiff er" and sup- 
presses number fluctuations relative to the ideal Bose 
gas [ 70ll84l . A complementary explanation is based on 
the observation that the condensate mode ardr is a low- 
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Figure 14: (Color online) First three moments of P{Nc) vs. 
total particle number (A*') of sGPe (open circles) vs. the- 
ory of Ref.|71| (filled green squares): (a) mean value (Nc) 
scaled to (A*') (condensate fraction), (b) relative standard de- 
viation a{Nc)/{Nc), (c) skewness or third centered moment, 
skew(Afc) = {{N^-{N^)f)/a^(N^). The total particle number 
N is calculated over the region \z\ < 2R. 



energy subsystem of the total field (represented by ip), 
where the non-condensate fraction can play the role of 
a particle reservoir. This suggests that the condensate 
subsystem can be described within a grand-canonical 
scheme even if the total atom number is fixed: for that it 
would be sufficient to consider a high enough tempera- 
ture so that a large number of non-condensate atoms is 
present. Indeed, the width cr(iVc) of the canonical count- 
ing statistics translates two physically different mech- 
anisms: on the one hand, the statistical uncertainty of 
the non-condensate (Bogoliubov) occupation numbers 
in the ncB expansion (exponentially distributed with 
mean N{Ek)), and on the other hand, the dynamical 
particle exchange with the non-condensate modes due 
to interactions, similar to what is done between system 
and bath within the sGPe. 

These considerations also suggest an explanation for 
the broader statistics that the sGPe method returns at 
low temperatures and small numbers [Fig|l3|d)]. It is 
symptomatic of the grand canonical ensemble which 
underlies the formulation of the sGPe, and leads to 
anomalously large number fluctuations for this nearly 
ideal gas. We have checked that the classical approxima- 
tion underlying the sGPe is not in error here: indeed, the 
counting statistics of the ncB method (canonical ensem- 
ble) is essentially the same when Bose-Einstein occupa- 
tion numbers are replaced by their classical (Rayleigh- 
Jeans) limit [Fig|l3|d), dotted blue curve]. 



Figure 14 shows the moments of the sGPe counting 
statistics as a function of the total particle number. This 
is compared to the theory of S&Co (canonical ensem- 
ble). We vary (N) over more than one order of mag- 
nitude, as a way to change the importance of particle 
interactions. The temperature is kept at a fixed ratio 
T/Tc = 0.23, where Tc = Tc{N) is the critical temper- 
ature for an ideal gas (see Eq.([lJ and Fig.|4|. The mean 
values [Fig[l4ja)] agree well between the ensembles, as 
expected [70 1, except perhaps at the smallest particle 
numbers. The standard deviation [Fig|T4|b)] is larger at 
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small (N) where one is closer to the ideal gas, but con- 
verges to S&Co theory for larger systems. In the third 
moment [Fig 14 c)], which measures the deviation from 
Gaussian statistics, we see that, at rather small numbers, 
the sGPe predicts a more symmetric distribution com- 
pared to the negatively skewed distribution obtained in 
the canonical ensemble. This suggests that in weakly 
interacting systems, the small non-condensate fraction 
cannot protect the condensate, like a "buffer", against 
the Gaussian noise in the stochastic d3mamics. The 
skewness builds up at higher particle numbers where 
also more modes are highly occupied, which is of course 
the regime where the sGPe should perform well. 



3. Fragmentation 

For T <^ Tcj,, one mode dominates the system, 
as many atoms are condensed and phase coherent 
(nearly pure condensate), whereas at higher temper- 
atures, many modes become appreciably occupied 
(quasi-condensate). This behaviour can be made quan- 
titative by considering the set of eigenvalues of the one- 
body density matrix {ip* {z)'4}{z')) , i.e., the average occu- 
pations Nk of the corresponding modes {k — 1, . . . , M). 
We recall that the eigenmodes of the one-body density 
matrix are distinct from those of the Hamiltonian, since 
the latter is not quadratic in the field. 

The modes for fc < 30 are shown for the system at 
T = 1.3T0 in the inset of Fig 15 The low-lying modes 
(fc < 10) share a significant fraction of the total occupa- 
tion; this is a consequence of the short-range phase co- 
herence or quasi-condensation in the system, and sim- 
ilar to a "fragmented" condensate where many modes 
share a macroscopic occupation [9 128J. A quantitative 
measure of how many modes contribute with a signif- 



icant occupation can be given in terms of the Schmidt 
number S defined as I.129J 
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Figure 15: (Color online) Schmidt number versus scaled tem- 
perature for the sGPe simulations. Inset: at the temperature 
T ^ 1.3 T^, the fractional occupation Nk/{N) for the 30 low- 
est modes with the largest occupation; condensate fraction (in 
= mode above) {Nc)/{N) ^ 0.2, and (A^) ^ 23 800. 
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S-' = Y.fl (51) 

fe=0 

where fk is the fractional mode occupation given by 



M 



(52) 



fc=0 



The Schmidt number versus temperature extracted from 
sGPe simulations is shown in the main plot of Figure 



15 In the limit of zero temperature, it tends towards 1: 
this is the signature of a pure Bose-Einstein condensate. 
For temperatures approaching T^, it increases quickly 
(a rather good fit is 5 w 1 + &{TlT^f, the black dotted 
line). At temperatures where the phase coherence length 
is shorter than the condensate size, one may expect a 
scaling S cx R{T)/L^{T) which would be slower than 
linear. Modes outside the condensate region therefore 
contribute as well. 



IV. CONDENSATION VS. QUASI-CONDENSATION 

Due to the ID nature of the system we consider, fluc- 
tuations in the density and phase are suppressed at 
different characteristic temperatures |3|. In this sec- 
tion, therefore, we highlight the distinction between 
the phase and density coherent portions of the gas, for 
which we will use the terminology "condensate" and 
"quasi-condensate", respectively. This is an important 
problem rn its own right, since stochastic theories (such 
as the sGPe UMllSll, or the pGPe [41 J) automatically gen- 
erate total densities of the field. The condensate mode is 
extracted from these either using bimodal fits (not com- 
monly done, but experimentally well-known), or via 
the Penrose-Onsager (PO) prescription. This leads to a 
'gap' between stochastic approaches and theories where 
a symmetry-breaking argument (or a variant of it) as- 
signs a special "condensate mode" from the outset [22 J. 

To investigate this issue further, we wish to estab- 
lish here a more direct link between the PO conden- 
sate and the quasi-condensate often calculated in low- 
dimensional systems; we do this by directly comparing 
the PO mode of the sGPe data to the ah initio prediction 
of the modified Popov theory of Refs.||6l[7l|92l outlined 
rn Sec |IIE| As we shall show, this link has the advan- 
tage of providing an approximate PO condensate den- 
sity without performing additional manipulations of the 
stochastic data like the diagonalization of the one-body 
density matrix. 
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A. Identifying the quasi-condensate 



In Fig 16 we compare the quasi-condensate calcu- 
lated from the modified Popov theory (within the lo- 
cal density approximation, see Sec II E I to the defini- 
tion | [53| below, based on the density correlation func- 
tion (Sec III C 2 1. In order to reveal the different physics 
contained in each of the approaches, in this figure, we 
have chosen a relatively high temperature (T w 1.3 ss 
0.63 Tc) than for the previous comparison (chemical po- 
tential and interaction constant are kept the same). The 
breakdown of the ncB initial state in that regime con- 
strains us to use only the sGPe stochastic data for this 
comparison. Due to the mapping between moments of 
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Figure 16: (Color online) Normalized density profiles show- 
ing: Total density (sGPe: Solid black, noisy; modified Popov: 
dashed, brown), quasi-condensate density (sGPe: solid or- 
ange, noisy; modified Popov: dashed, blue) and (phase coher- 
ent) condensate densities (sGPe PO: noisy, turquoise; modified 
Popov: dot-dashed /dotted, maroon). The modified Popov 
condensate density is shown dotted at small distances from 
the trap centre, where the relation n'^{z) = ndz) [see Eq.([56)] 
breaks down. The grey shaded region shows the Penrose- 
Onsager condensate density. The dashed red line shows n{z) — 
nqc{z) of the modified Popov theory. Here T — 1.3 = 
0.63 and (A^) ^ 23 800. There are no ncB data because the 
ncB expansion no longer works at this temperature. 

the Bose field operator to the stochastic field (Sec.|ll]l, one 
may extract a quasi-condensate density within these ap- 
proaches in the following way 

nliz) = 2(v3^(z)^(z))2 - {4'\z)i,{z)^^{z)i^iz)). (53) 



This definition has been put forward in Ref.flSO], and 
implemented in Refs. lll6[|122l . Its equivalent form 



,{z) = n{z)^2-gi^){z), 



(54) 



has also been used at lower temperatures with the 
aim of extracting the (conventional) condensate mode 
in a 3D system 1.1201 , where phase fluctuations were 



not expected to contribute significantly. The compar- 
ison of both the quasi-condensate (solid orange and 
dashed blue in Figjl6j and the total densities (solid 
black and dashed brown) between sGPe and modified 
Popov (respectively) gives very good agreement |64|. 
The quasi-condensate density profile is at this temper- 
ature clearly distinguishable from both the total den- 
sity and the PO condensate (greyed area). The physi- 
cal meaning of the quasi-condensate in modified Popov 
theory is thus that part of the system where density fluc- 
tuations are reduced such that (7(^^(2) « 1, as is typi- 
cal for a single-mode coherent state. The plateau with 
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IS, on 



5*^^^ (z) = 2 outside the quasi-condensate in Fig 
the other hand, characteristic for a "chaotic" (or multi- 
mode) field 175] [125] . Note that this definition of the 
quasi-condensate is immune to phase fluctuations by 
construction |3|. Figure 16 shows that the density cor- 
relations obtained within the sGPe [Eq.||54|] indeed cap- 
ture a quasi-condensate density consistent with that of 
Popov theory. 



B. Identifying the Penrose-Onsager condensate density 

The modified mean field theory of Andersen et al. 
f6 7 1 splits the system into a quasi-condensate and other 
modes and thus avoids the "problem" of assuming the 
existence of long-range phase coherence. This makes 
the theory valid in arbitrary dimensions and at all tem- 
peratures. The approach can also capture the (conven- 
tional, PO) condensate mode by calculating the long- 
range limit of the one-particle density matrix [92] (see 
also Ref. |112| ). In a homogeneous system, this leads to 
the definition 



Uc = lim nqc e' 



(55) 



This procedure recovers exactly the Popov results for 
quantum depletion in two and three dimensions, as 
pointed out in Ref .[92 1. 

We wish to adapt Eq.|[55) to the trapped case and con- 
struct the quantity 

= n(z)y^2-,g(2)(z)5W(0,z). (56) 

We expect that this agrees well with the PO condensate 
density ndz) for large \z\. As appears plausible on phys- 
ical grounds, Eq.|[56) involves a combination of density 
and phase correlation functions. This leads us to the sec- 
ond interesting feature of Figjl6] the density n'^{z) de- 
fined by Eq.l 56 1 (dot-dashed and dotted, maroon) coin- 
cides, at large \z\, with the condensate density nc{z) ob- 
tained from the PO analysis of the one-body density ma- 
trix (solid light blue and shaded). This is important, as 
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Figure 17: (Color onlme) Comparison between Eq. l [56) (dot-dashed, maroon) and the PO mode due to diagonalizing the density 
matrix (solid turquoise and shaded) at three temperatures. Also shown is the quasi-condensate density from Eq.fSi} (noisy orange 
curve), used to generate the dot-dashed maroon densities. The data here is extracted from the sGPe simulations. 



it illustrates that the PO procedure produces a conden- 
sate that coincides with that fraction of atoms for which 
both phase and density fluctuations are reduced, i.e. the 
'true' condensate. Equation | [56| l is also appealing from a 
practical point of view, as the extraction of the PO mode 
from stochastic data usually requires diagonalization of 
the one-body density matrix, which for very large sys- 
tems can become a significant computational task. In 
contrast, the correlation functions g^^^O,z) and 



that enter Eq.|[56) are very straightforward to calculate, 
making the analysis much quicker. 

The condensate density predicted by Eq.|[56| agrees 
well with the one obtained by Penrose-Onsager analy- 
sis also for a range of temperatures, as shown in Fig 17 
This figure also illustrates how the distinction between 
quasi-condensate and PO condensate becomes more no- 
ticeable as T T^. For all temperatures, Eq.|[56) gives a 
good approximation to the PO density, with an increas- 
ing difference evident only in the central region of the 
trap. The size of this region at T = 1.37^ corresponds 
roughly to the extent over which the non-condensate 



part of '{z) is positively correlated (see Eq.|43i and 



the dot-dashed red curve in Fig. |9|d)). The origin of the 
'spike', then, lies in the fact that the thermal component 
is coherent over this small region too, so it is not just the 
condensate which contributes to g^^^ (z) at short scales. 
Again from Fig. |9j it is interesting to note that the non- 
condensate contribution to 5*^^-' (z) is actually significant 
at larger \z\ where it reduces the condensate contribu- 
tion (dashed blue curve). 



C. Application: (quasi) condensate fraction 

To test the reliability of the condensate density sug- 
gested in Eq.|[56|, we have calculated the fraction of 
atoms in the (quasi-) condensate over a range of tem- 
peratures. The data are shown in Fig 18 The two upper 



curves give the quasi-condensate fraction which is sys- 
tematically larger. At low temperatures, T ^ T^, con- 
densate and quasi-condensate numbers are very simi- 
lar, as already noted in |44. ,116 , 131]. In this range, we 
also get a good agreement between stochastic data and 
the modified Popov theory. As T increases towards T^, 
however, the deviation between condensate and quasi- 
condensate increases. This is clearly due to the shrink- 
ing of the PO condensate at its borders, as seen in Fig 17 
The difference among the two (PO) condensate numbers 
is due to the central "spike" that follows from Eq.|[56|: 
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Figure 18: (Color online) (Id) Quasi-condensate and PO con- 
densate numbers from the sGPe (black circles) and modi- 
fied Popov theory (hollow green triangles). The sGPe results 
are based on the condensate mode given by PO analysis of 
the one-body density matrix (Sec 111B| and on Eq.|53} for the 



quasi-condensate. For the integration over riqciz), only the re- 
gion \z\ < R{T) is taken into accoimt. The modified Popov 
data are based on the quasi-condensate density described in 
Sec HE and on the condensate density n'^{z) given in Eq.l[56), 
where g^^\Q,z) is calculated from Eq.(|40). 
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Figure 19: (Color online) Top row (a, b): total atomic density for the sGPe (solid, black) at equilibrium and tWncB (a) before 
(dashed red) and (b) after (dashed blue) evolution via the GPe. The tWncB data contain the correction to normal order; the 
evolution period is 5000cj~^. Bottom row (c, d): condensate/ thermal density (sGPe:solid green/brown; tWncB: dashed), as 
obtained by Penrose-Onsager analysis [see Fig|7|, (c) before (dashed red) and (d) after (dashed blue) evolution of the ncB state 
over 5000cj~^. Upper curves: condensate density, lower curves with maxima at [z] < i?: thermal density. 



at z = 0, one has necessarily .g^^^O) = 1 and n^(0) = 
nqc{0), but diagonalizing the one-body density matrix 



of the sGPe data yields nc(0) < nqc(O) [Fig 17|. Despite 



this difference, the two condensate numbers agree rela- 
tively well and show the same qualitative temperature 
dependence. 

The approximate condensate density n^(z) of Eq.|56 
can be evaluated in a simpler way by using Eq.|40 
for the calculation of the phase correlation function 
(?^^-'(0, z). This means in practice that one only needs 
to extract the parameter R{T) from the modified Popov 
density profiles. We have already shown in Figj9|d) that 
this simplified expression works well to approximate 
5(i)(0,z). 



V. SLOW THERMALIZATION OF THE INITIAL STATE 

We have seen the importance of a consistent treatment 
of phase and density fluctuations in calculating accu- 
rately a finite-temperature initial state for our (quasi)- 
one-dimensional system. From the physical observables 
probed, the picture emerging so far is that equilibrium 
properties, like total densities, agree quite well over a 
range of temperatures well beneath and in the spa- 
tial range where highly occupied modes are dominant. 
The stochastic ensembles prepared by the two methods 
are considered in this section as an initial condition for 
the Gross-Pitaevskii equation, which will be used to de- 
scribe the subsequent dynamics. We take a tempera- 



ture w 0.48 where the ncB data is clearly not in equi- 
librium, and address the question whether these data 
evolve dynamically into a thermal equilibrium state. We 
shall find that even after a fairly long evolution time 
(> 5000 w^^), the system is not yet stationary. This may 
be related to the absence of thermalization in integrable 
homogeneous ID systems |[l08l[T09lll32]ri33l . 

The thermalization study presented here provides a 
link to other classical field methods. These approaches, 
for example the pGPe of Davis, Blakie and co-workers 
fiT], that employed by Berloff and Svistunov [54] and 
the approach of the Polish group [59 1, typically use a 
single field realization with suitably randomized initial 
conditions as an input to evolution under the GPe (with 
the possible addition of a projector [53 J), rather than an 
ensemble of initial states. The system then corresponds 
to a microcanonical ensemble, as both particle number 
and total energy are fixed in this scenario. The dynam- 
ics acquires an irreversible character by spatial or tem- 
poral coarse-graining [134]. This can be mapped to a 
Boltzmann equation that yields an irreversible evolution 
where the system thermalizes to an equilibrium state 
with Rayleigh-Jeans statistics ||4T||135I . 

Figure [19] shows the density profile and its resolution 
into condensate and thermal density, before (left panels) 
and after (right panels) GPe evolution. While the total 
densities (top row) agree well between the two meth- 
ods, the ncB data evolve towards a smaller condensate 
density (bottom row), suggesting the system to be ther- 
malized at a higher temperature. 
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Figure 20: (Color online) Left: cor\densate statistics at different evolution times under the GPe, at initial nominal temperature 
T « 0.48 Ttf,. Solid blue: tWncB data, dashed green: equilibrium S&co theory [71]. The sGPe data at equilibrium are not shown, 
as the condensate statistics essentially stays the same as in Fig[l3] Right: mean condensate number {Nc{t))for sGPe (upper, black 
line), tWncB (lower, red curve) at the temperatures of Table [ll} with temperature increasing from top (g) to bottom (i). 



The temporal evolution of the condensate statistics is 
illustrated in Fig.|20](left). This data was extracted as fol- 
lows: we obtain the one-body density matrix at the in- 
dicated times and get the condensate mode by Penrose- 
Onsager diagonalization. Projecting the ensemble of 
wave functions onto this mode, one gets a few snapshots 
of the evolving condensate statistics P{Nc, t). While this 
distribution shows essentially no variation for the sGPe 
initial data, the ncB case shows a significant evolution. 
The peak at (Nc) disintegrates quite rapidly, and the 
condensate is re-formed gradually, with a broader peak 
re-appearing from smaller to higher numbers. At the fi- 
nal time, the distribution does yet not appear to have 
reached a stable distribution. 



The right panel in Fig. 20 illustrates that over short 
time scales, the average condensate number oscillates 
for the ncB data, in contrast to the sGPe (note the shorter 
time scale compared to the left panels, in order to re- 
solve the oscillations in the ncB data). This oscillation 
at roughly 2 cj, whose amplitude increases with temper- 
ature, may be due to a nonlinear locking between the 
condensate mode and its low-lying, highly excited exci- 
tations. 

The evolution of the coherence functions is illustrated 
in Fig 21 for the 5^^^ and 5^^^ functions introduced in 



Sees III C 1 III C 2 The kinks in Fig. 21 a) could be due 
to a mode coupling between the condensate and low- 
energy excitations whose mode functions have nodes 
and are slightly broader. But it is unclear whether 
this picture may explain the strong density fluctuations 
(panel (d)). The average density is quite low for \z\ > R, 
which amplifies sampling errors, however. 



and |46l], we can use the correlation functions to mea- 
sure the temperature of the ensemble. This may not 
yield a consistent picture, since the system is not (yet) 
thermalized. At least, we can place bounds on the tem- 
perature range that the system might thermalize to, al- 
beit after some longer time. The results of this are sum- 
marized in Fig. 22 We probe the system properties at 
non-constant time intervals, in order to account for the 
possibility of periodic behaviour in the correlation func- 
tions. For example, the last two data points on the upper 
curve are separated by only lOoj"^ units of time. The 
lower dashed red line indicates the input temperature 
(0.48 T^), but both g'^' and g^-^' yield higher numbers. 
Notice that g'^^^ gives a higher temperature, as might be 
expected in our regime due to pronounced phase fluc- 
tuations. The temperature extracted from g^"^^ is more 
stable in time, which may suggest a faster damping 
rate for modes with density fluctuations. This would 
be consistent with Landau-Beliaev damping, see, e.g., 
Refs.r38l|59l. 



Recalling the discussion in Sees III CI III C 2 [Eqs.||42 



To summarize this discussion, we emphasize that in 
this example thermalization proceeds quite slowly. The 
ID character of the system that is nearly integrable prob- 
ably plays a role here, but also relevant are the large 
fluctuations that are present in the initial state produced 
by the ncB method. One may think of the 1/2 "quan- 
tum atoms" per mode that are included in the truncated 
Wigner sampling: the observed temperature increase is 
roughly comparable to the "classical thermalization" of 
these atoms, as simple estimates show [38 1. The differ- 
ent temperatures extracted from phase and density cor- 
relations, however, are probably related to the wrong ac- 
count of phase fluctuations: they are mis-interpreted in 
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Figure 21: (Color online) Top row: '{z) for the sGPe (solid black), (a) tWncB initially at {t = 0) (dastied red) and (b) tWncB 
after GPe evolution up to t = 5000tj^^ (dashed blue); bottom row: as per top row for g^^\z). In each plot, the vertical dashed 
line indicates R{T) at T = ASOIiuj. 



terms of non-condensate density, as illustrated in Fig 12 
This may be cured by formulating the stochastic scheme 
in terms of phase and density variables, instead of the 
non-condensate field ^±{z), similar to the analysis of 
Ref.L39J. 



VI. CONCLUSIONS 

We have analyzed the equilibrium properties of 
a weakly interacting, trapped quasi-one-dimensional 
Bose gas at finite temperatures, which we modelled as 
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Figure 22: (Color online) Temperature measurement using cor- 
relation functions as a function of GPe evolution time: (upper 
data) g^^\z), for < z < R/2 and Eq. (|42f (blue triangles); 
(lower data) g^^\0) and Eq. |46f (black squares). Initial tem- 
perature (dashed red); Best horizontal fit through p'^' (z) data 
(dotted brown) and through g'^'(O) data (dot-dashed green); 
temperatures corresponding to the horizontal lines are indi- 
cated. 



an effective one-dimensional system. The predictions 
of a number of independent finite-temperature theo- 
ries have been compared. We focussed in particular 
on two methods incorporating phase and density fluc- 
tuations in a stochastic manner: a number conserving 
Bogoliubov approach with stochastic sampling and the 
stochastic Gross-Pitaevskii equation. 

At low temperatures, we found average quantities, 
such as total density profiles, condensate fractions and 
first and second order spatial correlation functions to 
give good agreement between the two theories. These 
properties were additionally found to coincide with 
other theoretical predictions from the literature, which 
we used in order to 'benchmark' our findings. As 
higher temperatures were probed, though still within 
the regime T < T^, the ncB initial state was found to 
give predictions for equilibrium properties in disagree- 
ment with both the sGPe and the results of other, 'bench- 
mark' theories, including the modified Popov theory of 
[6], while the latter two showed good agreement. We 
attribute this failure to the use of the T = Bogoli- 
ubov spectrum in the ncB expansion, i.e. the conden- 
sate number is assumed equal to the total particle num- 
ber at all temperatures, and to the overestimation of the 
non-condensate density due to spurious contributions 
of phase fluctuations at higher orders. The 'point cloud' 
in Fig. 12 illustrates the enhanced phase fluctuations 
at higher temperatures, which mean that density fluc- 
tuations in modes above the condensate were no longer 
suppressed in the ncB method, as would be expected at 
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temperatures much below the 'degeneracy' temperature 
(Trf). A procedure taking condensate depletion into ac- 
count in a temperature dependent way, would likely im- 
prove the high-temperature behaviour of this method. 

We have also probed quantities involving higher- 
order statistical moments than just the density or its cor- 
relation, in particular the full distribution function of 
the condensate number At low temperatures, and for 
not too small atom numbers (where the assumption of 
a 'classical' occupation of modes fails), both ncB and 
sGPe were found to produce the correct statistics, in per- 
fect agreement to the theory of Svidzinsky and Scully 
IITTI . However, for small total particle numbers, and 
at low temperature, the sGPe results were found to be 
broader than the ncB statistics, whereas the latter were 
found to agree well with those of Ref . [71 1. The incorrect 
sGPe prediction at low particle numbers was attributed 
to the onset of anomalously large number fluctuations, 
familiar from the grand-canonical analysis of the ideal 
gas. As the importance of interactions within the system 
was increased, i.e. by increasing particle number with 
all other parameters fixed, the sGPe and Svidzinsky and 
Scully results were then found to match well. 

We further propagated the ncB data via the (ordinary) 
GPe, as done in the truncated Wigner approach. We 
found that this leads to the correct profiles for the total 
system density, but fails to predict all other features ac- 
curately, due to its attempt to thermalize to a higher tem- 
perature classical field. This thermalisation was found 
to take extremely long here, due to our Id system con- 
figuration. 

Finally, we have illustrated the conceptual differ- 
ence between the (phase-coherent) condensate and the 
(density-coherent) quasi-condensate. The former is usu- 
ally obtained by the Penrose-Onsager analysis of the 
one-body density matrix while the latter appears e.g. in 
the context of the modified Popov theory of Refs.j^lITI. 
Building on the identification of the (PO) condensate for 
a homogeneous system [92], we have provided an al- 
ternative, numerically very efficient means of extracting 
information about the condensate density that involves 
only first- and second-order correlation functions, as ob- 
tained from the sGPe simulations. Although some is- 
sues remain to be improved near the trap centre, con- 
densate density and fraction are perfectly matched to 
the conventional PO approach over a broad range of pa- 
rameters. We believe that this identification, along with 
the systematic benchmarking of observables to alterna- 
tive theories for finite-T Bose gases will provide a better 
understanding of the links between stochastic theories 
and thermodynamics based on mean field theories. 
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Appendix A: Sampling the number-conserving Bogoliubov 

state 

The condensate mode function (t)c{z) contains two 
terms 



<l>c{z) 



cj)^{z) + c^2{z)/N 

(l + ||02W/A^|P)l/2 



(Al) 



that arise in zero'th and second order of the expansion. 
The lowest-order contribution 0o(^) solves the station- 
ary GPe 



(A2) 



Note that the 'chemical potential' /i emerges here as the 
lowest eigenvalue of a nonlinear eigenproblem: it de- 
pends on the product gN and the trapping potential 
V{z). Equation | |A2| is conveniently solved by propa- 
gating the wave function in imaginary time. 

For the correction (t)2{z), one needs the non- 
condensate field %Ij^{z), and we return to it in Eq.ljASj. 
The non-condensate atoms populate Bogoliubov mode 
functions Uk{z) and Vk{z). These are eigenfunctions of 
the (projected) Bogoliubov-de Gennes operator 



Q* / Icq* 



c = 



i/Gp[2Af|'/'oP]-/i 



*2 



9N4>1 
-i?Gp[2iV|0o|']+A^ 



(A3) 



(A4) 



where Q (Q*) is the projector orthogonal to (t)o{z) (to 
(/)q(2:)), respectively. On the spatial grid, its matrix ele- 
ments are 



Az0nW0S(^')- 



(A5) 



The numerical diagonalization of Cq may be ap- 
proached simply with a Fourier grid method 11361, for 
example. We need in Eq.||6]l only those modes with 
Efe > that can be normalized to 

^2 X! ["^ {z)ui{z)~- vl {z)vi (z)] = Ski (A6) 
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where the sum represents the Bogoliubov scalar product 
on the spatial grid (spacing Az). For a grid of length L, 
the number of Bogoliubov modes is = L/Az — 1. 
The quantity A{{bk}) in Eq.||9} is calculated as 



-4({64) = E 



Az 

k,q 1 



E {Re(6fc 6* - 5kq(jl)vl{z)vq{z) 

z 

~Re{bk bq)ul{z)vq{z)} , 



(A7) 



where cr^. is the expectation value of |fefep, so that 
A{{bk}) averages to zero. This term encodes the 
Wigner correction for getting the (normally ordered) 
non-condensate particle number out of the semiclassi- 
cal Wigner functions, at the level of the corresponding 
number variances. 

The second-order correction 02(2;) to the condensate 
mode 00 (^) is due to thermal depletion. It contributes 
in the ncB expansion at the same order as the non- 
condensate modes to typical observables like the aver- 
age density and the one-body density matrix. It must 
be orthogonal to (j)o{z) (see ISOl ) and can hence be ex- 
panded over the Bogoliubov modes 



(A8) 



constructing an imaginary time evolution that leads to 
4)2{z) as a stationary solution. This is then plugged into 
Eq.| |Al| to complete the construction of the normalize 
condensate mode ipciz)- 

With this normalisation, the stochastic matter wave 
field V'C^) gives access to the total particle density n{z) 
as [cf. Eq.p] 



The coefficients are found by solving an inhomoge- 
neous linear equation for the Bogoliubov-de Gennes op- 
erator. For the convenience of the reader, we reproduce 
here Eqs.(71)-(75) provided in Ref.|46|. The Wigner field 
4>± represents a non-condensate field operator denoted 
A in Ref. ll46l . The second-order correction ^2 solves the 
stationary equation 

Q(ffGp[2iV|0o|'] - A')02 + gNQcj,lr2 = -QR (A9) 
with a "source term" 

R{x) = -giV|0o(x)|'0o(x)(l + (iVth» (AlO) 
+2gN{k\x)k{x))Mx) + gN<l,l{x){k{x)k{x)) 

-gN [ dy \My)\'{i^Hy)My) + ro(.y)Hy))m) 



where iVth = Jda; A^^ {x)k{x) is the operator for the non- 
condensate atom number. The expectation values of A 
are translated into Wigner averages of tjj± according to 
the symmetrization rule 

{k\x)k{x')) - (V.1(x)V'^(x'))m/ - ^Q(x,x') (All) 

where the projector Q(x, a;')appears because the fields 
live in the subspace orthogonal to the zero'th order 
condensate mode (t)Qiz). Equation |A9i is solved by 



n-{z) = {\ip{z)\'^)w - nq 



{All) 



where is the quantum density It is hence normal- 
ized such that {\\ip\\^)w = N + M/2. The statistics of 
the condensate atom number Nc is determined by the 
non-condensate field V'i/ via Eq.||9|. 

The propagation in time of the matter field ip{z, t) is 
found by solving the GPe l|2j. This turns out to be more 
accurate, in particular at long times, than to propagate 
separately the condensate mode and using the time- 
dependent versions of the equations they solve in the re- 
spective orders of the expansion in {5N/Ny^^ (see, e.g., 
113 7|| ). All numerical simulations in this paper (for both 
methods) are based on the Crank-Nicholson method for 
time stepping with the results averaged over at least 
1000 realisations of the initial conditions. 



Appendix B: Canonical counting statistics of a Bose gas 

The group of Scully and co-workers (S&Co) has devel- 
oped theoretical models to calculate the counting statis- 
tics P(A*'c) of a Bose condensate, building on the canoni- 
cal ensemble where the operators Nc and Nth must sum 
up to the (fixed) total number N. Therefore, the count- 
ing statistics is the "mirror image" of the probability dis- 
tribution P(iVth). The latter can be calculated when the 
non-condensate number iVth splits into a sum of sta- 
tistically independent terms. In an ideal Bose gas, this 
would be the occupation numbers fik of single-particle 
modes (quantum number k). For a weakly interacting, 
homogeneous Bose gas, the non-condensate atoms are 
clearly those in non-zero momentum modes, and hence 

p>0 

+ 2upVpiblblp + bpb_p) + 2vl} (Bl) 

where Up and Vp are real-valued Bogoliubov coefficients 
[Eq.(258) of Ref.|70|], normalized to - = 1, and 
where the bp are the bosonic operators for Bogoliubov 
quasi-particles. This Bogoliubov spectrum depends on 
the condensate occupation number Nc- S&Co make the 
approximation that the dependence is weak enough so 
that Nc can be replaced by its average value {Nc) (i.e., by 
the first moment of P{Nc)). This is consistent if P{Nc) 
is narrow enough. The operator identity ||B1| is also 
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based on the assumption of a negligibly small probabil- 
ity P{Nc ~ 0) of finding the condensate mode empty 
[Eq.(210) of Ref.EQl]. 

The quasi-particle operators are constructed such that 
the sum | |B1| contains mutually commuting operators 
(only momenta p and — p are correlated), and the prob- 
ability distribution P(iVth) can be found with standard 
techniques. The following moments are found, for ex- 
ample 

(iVth) = ^{TVpU^ + (iVp + (B2) 
^^'(^^th) - J2^Np{Np + l)[l + 8ulvl] + 2ulvl} (B3) 

X [1 + Wulvl] + {2Np + l)4ulvl} (B4) 

where Np = N{ep) is the Bose-Ernstern statistics and 
where ^3 = ((A^th — (-^th))^) is the third central mo- 
ment. Its nonzero value is a clear indication of non- 
Gaussian statistics. These results are obtained within the 
Bogoliubov approximation with a weak thermal frac- 
tion, (A^th) ^ N, hence at low temperatures. 

Svidzinsky and Scully \71\ have generalized this ap- 
proach to any temperature, using a rate equation Ansatz 
similar to quantum laser theory |74|. The growth and 
loss rates for the condensate depend on Nc and are de- 
scribed by a rational (Fade) approximation. This leads 
to a closed formula for the stationary P{Nc). The Fade 
parameters are matched to the low-temperature limit 

-|B4l. 



and can be expressed in terms of the moments I B2 - 
The equation (Nc) = N — {Nth) must be solved iter- 
atively since the Bogoliubov modes in Eq.(jB2j depend 
themselves on (Nc). It has been pointed out in Ref.[72| 
that the dependence of the Bogoliubov spectrum on the 
condensate number, ep{Nc) ^ ep{{Nc)), actually leads 
to an observable difference in the counting statistics, al- 
though it makes the calculations much more involved. 

We calculate the counting statistics in Figs 13 14 us- 
ing the theory of Ref.|71 1 and making the identification 
Up I— > ||w/c|P to the Bogoliubov modes in a trap. The mo- 



B4l 



ments of P{Nc) in Fig 14 are calculated from Eqs.|B2- 
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